{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys,os,inspect\n",
    "if hasattr(sys.modules[__name__], '__file__'):\n",
    "    _file_name = __file__\n",
    "else:\n",
    "    _file_name = inspect.getfile(inspect.currentframe())\n",
    "CURRENT_FILE_PATH = os.path.dirname(_file_name)\n",
    "sys.path.append(os.getcwd()+\"/../neuronVis\")\n",
    "sys.path.append(os.getcwd()+\"/../resource\")\n",
    "from cx_tool import *\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Information generation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_path = '../resource/project_insula/'\n",
    "create_dir(data_path)\n",
    "_flatmap = cv.imread(data_path+'flatmapedge2.tif',1)\n",
    "cortex_region = ['AI','GU','VISC']\n",
    "cortex_region_to_subregion,cortex_subregion_to_region = get_subregions(cortex_region)\n",
    "target_region_hierarchy = {\n",
    "        'Cortex':['FRP','PL','ILA','ORB','ACA','AI','GU','VISC','PTLp','RSP','TEa','ECT','PERI','MOs','MOp','SSp','SSs','AUD','VIS'],\n",
    "        'OLF':['MOB','AOB','AON','TT','DP','PIR','NLOT','COA','PAA','TR'],\n",
    "        'HPF':['HIP','RHP'],\n",
    "        'Cortex-sp':['CTXsp'],\n",
    "        'CNU':['STRd','STRv','LSX','sAMY','PAL'],\n",
    "        'TH':['VENT','SPF','SPA','PP','GENd','LAT','ATN','MED','MTN','ILM','RT','GENv','EPI'],\n",
    "        'HY':['PVZ','PVR','MEZ','LZ','ME'],\n",
    "        'MB':['VTA','SNr','SCs','SCm','NB','SAG','PBG','IC','RR','MRN','RN','PAG','PRT','CUN','PPN','RAmb','III','MA3','EW','VTN'],\n",
    "        'P':['PSV','PB','PCG','PG','PRNc','SUT','TRN','V','P5','CS','LDT','NI','PRNr','RPO','SOC','NLL','DTN','PDTg','SG','I5','SLC','SLD'],\n",
    "        'MY':['CN','DCN','ECU','NTB','NTS','SPVC','SPVI','SPVO','Pa5','VII','GRN','IO','IRN','LIN','LRN','MARN',\n",
    "        'MDRNd','MDRNv','PARN','PGRNd','PGRNl','PHY','VNC','VI','AMB','DMX','ICB','PAS','PRP','PPY','x','XII','RM','RPA','RO'],\n",
    "        'CB':['CB'],\n",
    "        'FT':['py','pyd']\n",
    "    }\n",
    "target_region = []\n",
    "for i in target_region_hierarchy:\n",
    "    for j in target_region_hierarchy[i]:\n",
    "        target_region.append(j)\n",
    "if not os.path.exists(data_path+'cortex_neuron.npy'):\n",
    "    cortex_neuron = {}\n",
    "    for i in cortex_region:\n",
    "        for j in iondata.getNeuronListBySomaBrainRegion(i):\n",
    "            if len(j['sampleid']+j['name'])!=13:\n",
    "                continue\n",
    "            cortex_neuron[j['sampleid']+j['name']] = {}\n",
    "    if not os.path.exists(data_path+'ITCTPT_mask.npy'):\n",
    "        mask = np.where(get_region_mask('CTX')>0,1,0)\n",
    "        mask[get_region_mask('CNU').astype(bool)] = 1\n",
    "        mask[get_region_mask('TH').astype(bool)] = 2\n",
    "        mask[get_region_mask('HY').astype(bool)] = 3\n",
    "        mask[get_region_mask('MB').astype(bool)] = 3\n",
    "        mask[get_region_mask('HB').astype(bool)] = 3\n",
    "        mask[get_region_mask('CB').astype(bool)] = 3\n",
    "        np.save(data_path+'ITCTPT_mask.npy',mask)\n",
    "    else:\n",
    "        mask = np.load(data_path+'ITCTPT_mask.npy')\n",
    "    def get_nearest_point(p):\n",
    "        p = np.array(p)\n",
    "        if mask[p[0],p[1],p[2]]!=0:\n",
    "            return mask[p[0],p[1],p[2]]\n",
    "        else:\n",
    "            r = 2\n",
    "            while 1:\n",
    "                start = p-r\n",
    "                end = p+r+1\n",
    "                start = np.where(start<0,0,start)\n",
    "                tmp = Counter(mask[start[0]:end[0],start[1]:end[1],start[2]:end[2]].flatten().tolist()).most_common(2)\n",
    "                if len(tmp)>1:\n",
    "                    if tmp[0][0]==0:\n",
    "                        return tmp[1][0]\n",
    "                    else:\n",
    "                        return tmp[0][0]\n",
    "                r += 5\n",
    "    def f(neuron,mode=0):\n",
    "        neuron_tree = iondata.getNeuronTreeByID(neuron['sampleid'],neuron['name'])\n",
    "        info = {'A':[],'B':[],'C':[],'D':[]}\n",
    "        info2 = {'A':[],'B':[],'C':[],'D':[]}\n",
    "        for terminal in neuron_tree.terminals:\n",
    "            if terminal.type==3:\n",
    "                continue\n",
    "            tmp = (np.array(terminal.xyz)/10).astype(int)\n",
    "            tmp = [np.clip(tmp[i],0,mask.shape[i]-1) for i in range(3)]\n",
    "            if mode:\n",
    "                region = get_nearest_point(tmp)\n",
    "            else:\n",
    "                region = mask[tmp[0],tmp[1],tmp[2]]\n",
    "            if region==1:\n",
    "                info['A'].append(terminal.index)\n",
    "                info2['A'].append(terminal.xyz)\n",
    "            elif region==2:\n",
    "                info['B'].append(terminal.index)\n",
    "                info2['B'].append(terminal.xyz)\n",
    "            elif region==3:\n",
    "                info['C'].append(terminal.index)\n",
    "                info2['C'].append(terminal.xyz)\n",
    "            else:\n",
    "                info['D'].append(terminal.index)\n",
    "                info2['D'].append(terminal.xyz)\n",
    "        cls = 'Others'\n",
    "        a = [len(info[i]) for i in info.keys()]\n",
    "        if a[0]!=0 and a[1]==0 and a[2]==0:\n",
    "            cls = 'IT'\n",
    "        if a[1]!=0 and a[2]==0:\n",
    "            cls = 'CT'\n",
    "        if a[2]!=0:\n",
    "            cls = 'PT'\n",
    "        return [cls,info,info2]\n",
    "    for i in tool_bar(cortex_neuron,desc='iteration on neurons (generating information) ...',leave=False,colour='#55AF7B'):\n",
    "        control_console_output(0)\n",
    "        neuron_property = iondata.getNeuronPropertyByID(i[:-7],i[-7:])\n",
    "        cortex_neuron[i]['soma_position'] = neuron_property['somapoint']\n",
    "        cortex_neuron[i]['soma_hemisphere'] = 'left' if neuron_property['somapoint'][2]<5700 else 'right'\n",
    "        a,b = region_layer_seperate(neuron_property['somaregion'])\n",
    "        cortex_neuron[i]['cortex_region'] = cortex_subregion_to_region[neuron_property['somaregion']]\n",
    "        cortex_neuron[i]['cortex_subregion'] = a[len(cortex_neuron[i]['cortex_region']):]\n",
    "        cortex_neuron[i]['cortex_layer'] = b\n",
    "        j = compute_flatten_position_and_depth([get_mirrored_point(cortex_neuron[i]['soma_position'])])[0]\n",
    "        cortex_neuron[i]['soma_depth'] = j[0]\n",
    "        cortex_neuron[i]['soma_outter_position'] = j[1]\n",
    "        cortex_neuron[i]['soma_inner_position'] = j[2]\n",
    "        cortex_neuron[i]['soma_flatten_position'] = j[3]\n",
    "        cortex_neuron[i]['type'] = f({'sampleid':i[:-7],'name':i[-7:]},mode=1)[0]\n",
    "        control_console_output(1)\n",
    "    for i in tool_bar(cortex_neuron,desc='generating neuron information...',leave=False):\n",
    "        cortex_neuron[i]['axonal_projection_length_combined'] = get_axonal_statics_by_regions(i,target_region,'projectregion')\n",
    "        cortex_neuron[i]['axonal_projection_length_left'] = get_axonal_statics_by_regions(i,target_region,'projectleftregion')\n",
    "        cortex_neuron[i]['axonal_projection_length_right'] = get_axonal_statics_by_regions(i,target_region,'projectrightregion')\n",
    "        cortex_neuron[i]['axonal_projection_terminal_combined'] = get_axonal_statics_by_regions(i,target_region,'terminalregion')\n",
    "        cortex_neuron[i]['axonal_projection_terminal_left'] = get_axonal_statics_by_regions(i,target_region,'terminalleftregion')\n",
    "        cortex_neuron[i]['axonal_projection_terminal_right'] = get_axonal_statics_by_regions(i,target_region,'terminalrightregion')\n",
    "    np.save(data_path+'cortex_cluster_object.npy',cortex_neuron)\n",
    "else:\n",
    "    cortex_neuron = np.load(data_path+'cortex_neuron.npy',allow_pickle=True).item()\n",
    "cortex_neuron_itctpt = {'IT':[],'CT':[],'PT':[],'Others':[]}\n",
    "for i in cortex_neuron:\n",
    "    cortex_neuron_itctpt[cortex_neuron[i]['type']].append(i)\n",
    "len_all = len(cortex_neuron.keys())\n",
    "len_it = len(cortex_neuron_itctpt['IT'])\n",
    "len_ct = len(cortex_neuron_itctpt['CT'])\n",
    "len_pt = len(cortex_neuron_itctpt['PT'])\n",
    "len_ot = len(cortex_neuron_itctpt['Others'])\n",
    "print('Currently there are',len_all,'neurons in total.\\n',\n",
    "    'The IT neurons are',len_it,' (',len_it/len_all*100,'%).\\n',\n",
    "    'The CT neurons are',len_ct,' (',len_ct/len_all*100,'%).\\n',\n",
    "    'The PT neurons are',len_pt,' (',len_pt/len_all*100,'%).\\n',\n",
    "    'The Other neurons are',len_ot,' (',len_ot/len_all*100,'%).')\n",
    "fontsize_ticks = 20\n",
    "fontsize_label = 25\n",
    "fontsize_title = 30"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Region transformation (A/M/P version)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Xiao's cortex region transformation\n",
    "IC_a_m = -1000*1.21+5472\n",
    "IC_m_p = -1000*0.01+5472\n",
    "insular_neuron_name = []\n",
    "for i in cortex_neuron:\n",
    "    insular_neuron_name.append(i)\n",
    "region_transform_region = ['AIa','AIm','AIp','GU','VISC']\n",
    "region_transform_color = {'AIa':'#DB4537','AIm':'#55AF7B','AIp':'#4286F3','GU':'#95dfe3','VISC':'#FAC230'}\n",
    "region_transform = {}\n",
    "for i in insular_neuron_name:\n",
    "    if cortex_neuron[i]['cortex_region']=='AI':\n",
    "        if cortex_neuron[i]['soma_position'][0]<IC_a_m:\n",
    "            region = 'AIa'\n",
    "        else:\n",
    "            if cortex_neuron[i]['soma_position'][0]<IC_m_p:\n",
    "                region = 'AIm'\n",
    "            else:\n",
    "                region = 'AIp'\n",
    "        region_transform[i] = region\n",
    "    else:\n",
    "        region_transform[i] = cortex_neuron[i]['cortex_region']\n",
    "flatmap = cv.imread('../resource/flatmapedge2.tif',1)\n",
    "fig = plt.figure(figsize=(5,5))\n",
    "ax = fig.add_subplot(1,1,1)\n",
    "group = {}\n",
    "for i in insular_neuron_name:\n",
    "    region = region_transform[i]\n",
    "    position = cortex_neuron[i]['soma_flatten_position']\n",
    "    if region not in group:\n",
    "        group[region] = []\n",
    "    group[region].append(position)\n",
    "ax.matshow(flatmap)\n",
    "for i in group:\n",
    "    tmp = np.array(group[i])*2\n",
    "    ax.scatter(tmp[:,0],tmp[:,1],marker='o',color=region_transform_color[i],s=0.1)\n",
    "    ax.scatter([],[],marker='o',color=region_transform_color[i],s=20,label=i) \n",
    "ax.set_xticks([],[])\n",
    "ax.set_yticks([],[])\n",
    "#ax.set_axis_off()\n",
    "ax.set_xlabel('Lat. to Med. axis',loc='center',fontsize=fontsize_label*0.6)\n",
    "ax.set_ylabel('Ant. to Post. axis',rotation=-90,loc='center',fontsize=fontsize_label*0.6).set_verticalalignment('top')\n",
    "ax.set_title('Soma distribution of the insular neurons',fontsize=fontsize_title*0.5)\n",
    "ax.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = []\n",
    "y = []\n",
    "for i in insular_neuron_name:\n",
    "    x.append(get_mirrored_point(cortex_neuron[i]['soma_position']))\n",
    "    if 1:\n",
    "        if cortex_neuron[i]['soma_position'][0]<IC_a_m:\n",
    "            y.append([0.859,0.267,0.216])\n",
    "        else:\n",
    "            if cortex_neuron[i]['soma_position'][0]<IC_m_p:\n",
    "                y.append([0.059,0.616,0.345])\n",
    "            else:\n",
    "                y.append([0.258,0.52,0.957])\n",
    "neuronvis = nv.neuronVis()\n",
    "neuronvis.render.setBackgroundColor((1.0,1.0,1.0,1.0))\n",
    "draw_points(neuronvis,x,y)\n",
    "neuronvis.addRegion('MO',[1,0,0])\n",
    "neuronvis.addRegion('SSs',[0,1,0])\n",
    "#neuronvis.addRegion('VISC',[0,0,1])\n",
    "neuronvis.render.setPointSize(5)\n",
    "neuronvis.render.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Anterior to Striatum dorsal and nucleus accumbens"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mask = get_region_mask('STRd')\n",
    "mask = np.where(get_region_mask('ACB')==1,2,mask)\n",
    "mask = np.where(get_region_mask('GPe')==1,3,mask)\n",
    "background,background_header = nrrd.read('background.nrrd')\n",
    "background = background.transpose((2,1,0))\n",
    "tmp = background*1\n",
    "del background\n",
    "background = tmp\n",
    "projection_mask = mask*0\n",
    "external_point = []\n",
    "c = 0\n",
    "for neuron in tool_bar(insular_neuron_name):\n",
    "    if region_transform[neuron]!='AIa':\n",
    "        continue\n",
    "    c += 1\n",
    "    neuron_tree = iondata.getNeuronTreeByID(neuron[:-7],neuron[-7:])\n",
    "    for point in neuron_tree.points:\n",
    "        if point.type==3:\n",
    "            continue\n",
    "        x = np.round(np.array(point.xyz)/10).astype(int)\n",
    "        try:\n",
    "            projection_mask[x[0],x[1],x[2]] += 1\n",
    "        except:\n",
    "            external_point.append(x)\n",
    "import numpy as np\n",
    "projection_mask = np.load('projection_mask.npy')\n",
    "def subcube(x,center,r):\n",
    "    slc = tuple(slice(max(c - r, 0), c + r + 1) for c in center)\n",
    "    subcube = x[slc]\n",
    "    padded_subcube = np.pad(subcube,[(max(r-center[i],0),max(center[i]+r-j+1,0)) for i,j in enumerate(x.shape)],mode='constant')\n",
    "    return padded_subcube\n",
    "\n",
    "\n",
    "size = 5\n",
    "sigma = 2.0\n",
    "\n",
    "filter_3d = np.fromfunction(\n",
    "    lambda x, y, z: (1/(2*np.pi*sigma**2)) * np.exp(-((x-(size-1)//2)**2 + (y-(size-1)//2)**2 + (z-(size-1)//2)**2) / (2*sigma**2)),\n",
    "    (size, size, size), dtype=float\n",
    ")\n",
    "\n",
    "# Normalize the filter\n",
    "filter_3d /= np.sum(filter_3d)\n",
    "#import tensorflow as tf\n",
    "#a = projection_mask.astype('float16').reshape(1,projection_mask.shape[0],projection_mask.shape[1],projection_mask.shape[2],1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "projection_mask_shrink = projection_mask.reshape(projection_mask.shape[0]//2,2,projection_mask.shape[1]//2,2,projection_mask.shape[2]//2,2)\n",
    "projection_mask_shrink = np.mean(projection_mask_shrink,axis=(1,3,5))\n",
    "projection_mask2 = projection_mask_shrink*0\n",
    "for i in tool_bar(range(projection_mask_shrink.shape[0])):\n",
    "    for j in range(projection_mask_shrink.shape[1]):\n",
    "        for k in range(projection_mask_shrink.shape[2]):\n",
    "            projection_mask2[i,j,k] = np.sum(subcube(projection_mask_shrink,[i,j,k],2)*filter_3d)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.ndimage import gaussian_filter\n",
    "projection_mask_processed = gaussian_filter(projection_mask*10,sigma=2)\n",
    "AIa_CNU_path = data_path+'CNU/'\n",
    "create_dir(AIa_CNU_path)\n",
    "data_frames = np.log(projection_mask_processed+1)\n",
    "try:\n",
    "    data_frames /= np.max(data_frames)\n",
    "except:\n",
    "    pass\n",
    "selected_region_frames = mask/5\n",
    "background_frames = np.log(background+1)**2.5\n",
    "try:\n",
    "    background_frames /= np.max(background_frames)\n",
    "except:\n",
    "    pass\n",
    "for index in range(150,850,25):\n",
    "    bregma = np.round((5472-index*10)/1000,2)\n",
    "    data_frame = data_frames[index]\n",
    "    selected_region_frame = selected_region_frames[index]\n",
    "    background_frame = background_frames[index]\n",
    "    data_mask = np.where(data_frame>0,0,1)\n",
    "    selected_region_mask = np.where(selected_region_frame>0,0,1)\n",
    "    merged_frame = np.array([data_frame,data_mask*selected_region_frame*1,data_mask*background_frame*1]).transpose(1,2,0)\n",
    "    plt.matshow(merged_frame)\n",
    "    plt.title('Bregma = '+str(bregma),fontsize=15)\n",
    "    plt.axis('off')\n",
    "    plt.savefig(AIa_CNU_path+str(bregma)+'.svg',format='svg')\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Allen version analysis (all)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tmp = []\n",
    "for i in cortex_neuron:\n",
    "    tmp.append(i)\n",
    "insular_neuron_name = tmp\n",
    "cluster_path = data_path+'cluster_allen/'\n",
    "create_dir(cluster_path)\n",
    "if_combined = False\n",
    "if if_combined:\n",
    "    neuron_axonal_length_by_region_log = []\n",
    "    for i in insular_neuron_name:\n",
    "        neuron_axonal_length_by_region_log.append([cortex_neuron[i]['axonal_projection_length_combined']])\n",
    "    neuron_axonal_length_by_region_log = np.array(neuron_axonal_length_by_region_log).reshape(-1,len(target_region))\n",
    "    for i in target_region_hierarchy['FT']:\n",
    "        neuron_axonal_length_by_region_log[:,target_region.index(i)] *= 2\n",
    "else:\n",
    "    neuron_axonal_length_by_region_log = []\n",
    "    for i in insular_neuron_name:\n",
    "        if cortex_neuron[i]['soma_hemisphere']=='left':\n",
    "            neuron_axonal_length_by_region_log.append([cortex_neuron[i]['axonal_projection_length_left'],cortex_neuron[i]['axonal_projection_length_right']])\n",
    "        else:\n",
    "            neuron_axonal_length_by_region_log.append([cortex_neuron[i]['axonal_projection_length_right'],cortex_neuron[i]['axonal_projection_length_left']])\n",
    "    neuron_axonal_length_by_region_log = np.array(neuron_axonal_length_by_region_log).reshape(-1,2*len(target_region))\n",
    "    for i in target_region_hierarchy['FT']:\n",
    "        neuron_axonal_length_by_region_log[:,target_region.index(i)] *= 2\n",
    "        neuron_axonal_length_by_region_log[:,target_region.index(i)+len(target_region)] *= 2\n",
    "    neuron_axonal_length_by_region_log = np.log(neuron_axonal_length_by_region_log+1)\n",
    "if not os.path.exists(cluster_path+'cluster_color_20.npy'):\n",
    "    cluster_color = np.array(distinctipy.get_colors(20,background_colors.tolist()))\n",
    "    np.save(cluster_path+'cluster_color_20.npy',cluster_color)\n",
    "else:\n",
    "    cluster_color = np.load(cluster_path+'cluster_color_20.npy')\n",
    "cluster_color_hex = []\n",
    "for i in cluster_color:\n",
    "    cluster_color_hex.append(rgb_to_hex(np.round(i*255).astype(int)))\n",
    "print('The colors of clustering:')\n",
    "distinctipy.color_swatch(cluster_color)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hierarchy.set_link_color_palette(cluster_color_hex)\n",
    "color_threhold = 300\n",
    "fig = plt.figure(figsize=(10,5))\n",
    "if not os.path.exists(cluster_path+'cluster_object2.npy'):\n",
    "    cluster_object = hierarchy.dendrogram(hierarchy.linkage(neuron_axonal_length_by_region_log, 'ward'),color_threshold=color_threhold,distance_sort=True)\n",
    "    np.save(cluster_path+'cluster_object.npy',cluster_object)\n",
    "else:\n",
    "    cluster_object = np.load(cluster_path+'cluster_object.npy',allow_pickle=True).item()\n",
    "plt.axhline(color_threhold, linestyle='dashed', color='#AAAAAA')\n",
    "plt.title(\"Hierachy of the clustering\",fontsize=20)\n",
    "plt.locator_params(axis='y', nbins=3)\n",
    "plt.xticks([],[])\n",
    "plt.xlabel('neurons = '+str(len(insular_neuron_name)),fontsize=15)\n",
    "plt.ylabel('Cluster distance',fontsize=15)\n",
    "plt.savefig(cluster_path+'cluster_hierachy.svg',format='svg')\n",
    "plt.show()\n",
    "cluster_object_parse = [[],[]]\n",
    "neuron_to_cluster = {}\n",
    "for i,j in zip(cluster_object['leaves'],cluster_object['leaves_color_list']):\n",
    "    neuron = insular_neuron_name[i]\n",
    "    projection_matrix = neuron_axonal_length_by_region_log[i]\n",
    "    cluster_object_parse[0].append(cluster_color_hex.index(j))\n",
    "    cluster_object_parse[1].append(projection_matrix)\n",
    "    neuron_to_cluster[neuron] = cluster_color_hex.index(j)\n",
    "cluster_object_parse[1] = np.array(cluster_object_parse[1])\n",
    "cluster_object_parse_visualization = np.zeros(shape=(len(insular_neuron_name),1+len(target_region),3))\n",
    "cluster_object_parse_visualization[:,0,:] = np.array([cluster_color[i] for i in cluster_object_parse[0]])\n",
    "tmp = np.array([cluster_object_parse[1][:,:len(target_region)],cluster_object_parse[1][:,:len(target_region)]*0,cluster_object_parse[1][:,len(target_region):]]).transpose((1,2,0))\n",
    "tmp /= np.max(tmp)\n",
    "cluster_object_parse_visualization[:,1:,:] = tmp*1\n",
    "fig = plt.figure(figsize=(20,60))\n",
    "ax = fig.add_subplot()\n",
    "ax.matshow(cluster_object_parse_visualization.transpose((1,0,2)))\n",
    "_ = ax.set_xticks([],[])\n",
    "ax.set_xlabel('neurons = '+str(len(insular_neuron_name)),fontsize=fontsize_label)\n",
    "cluster_object_parse_visualization_yticks = [[0],['class']]\n",
    "cluster_object_parse_visualization_yticks[0] += [1+i for i in range(len(target_region))]\n",
    "cluster_object_parse_visualization_yticks[1] += target_region\n",
    "_ = ax.set_yticks(cluster_object_parse_visualization_yticks[0],cluster_object_parse_visualization_yticks[1],fontsize=fontsize_ticks)\n",
    "plt.savefig(cluster_path+'cluster_projection_matrix.svg',format='svg')\n",
    "plt.show()\n",
    "tsne = TSNE(n_components=2, perplexity=35, random_state=0,n_iter=10000)\n",
    "X_tsne = tsne.fit_transform(neuron_axonal_length_by_region_log)\n",
    "tmp = {}\n",
    "for i,j in enumerate(insular_neuron_name):\n",
    "    _cluster = cluster_color_hex.index(cluster_object['leaves_color_list'][cluster_object['leaves'].index(i)])\n",
    "    if _cluster not in tmp:\n",
    "        tmp[_cluster] = []\n",
    "    tmp[_cluster].append([j,cortex_neuron[j]['soma_flatten_position'],X_tsne[i]])\n",
    "_flatmap = cv.imread('data/flatmapedge2.tif',1)\n",
    "fig = plt.figure(figsize=(10,10))\n",
    "flatmap = _flatmap*1\n",
    "for i in range(len(list(tmp.keys()))):\n",
    "    tmp2 = []\n",
    "    tmp3 = []\n",
    "    for j in tmp[i]:\n",
    "        if np.sum(abs(np.array(j[1])))!=0:\n",
    "            cv.circle(flatmap,[int(j[1][0])*2, int(j[1][1])*2],5,cluster_color[i]*255,-1)\n",
    "        tmp2.append(j[2])\n",
    "        tmp3.append(j[0])\n",
    "    tmp2 = np.array(tmp2)\n",
    "    plt.scatter(tmp2[:,0], tmp2[:,1],c=cluster_color_hex[i],s=30,marker='o',label=str(i+1))\n",
    "    save_scene(tmp3,cluster_color[i],cluster_path+'cluster_'+str(i+1))\n",
    "plt.legend()\n",
    "plt.title(\"t-SNE Visualization of the clusters\",fontsize=fontsize_title)\n",
    "plt.axis('off')\n",
    "plt.savefig(cluster_path+'cluster_t-SNE.svg',format='svg')\n",
    "plt.show()\n",
    "fig = plt.figure(figsize=(10,10))\n",
    "plt.imshow(flatmap)\n",
    "plt.xlim(300,2700)\n",
    "plt.ylim(2700,300)\n",
    "plt.title(\"The soma distribution of the clusters\",fontsize=fontsize_title)\n",
    "plt.axis('off')\n",
    "plt.savefig(cluster_path+'cluster_distribution.svg',format='svg')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(15,5))\n",
    "number_clusters = []\n",
    "for i in neuron_to_cluster:\n",
    "    if neuron_to_cluster[i] not in number_clusters:\n",
    "        number_clusters.append(neuron_to_cluster[i])\n",
    "number_clusters = len(number_clusters)\n",
    "kernel_size = 15\n",
    "kernel_sigma = 35\n",
    "for cluster in range(number_clusters):\n",
    "    distribution = np.zeros((1800,))\n",
    "    distribution_index = np.array([i for i in range(1800)])\n",
    "    for i in insular_neuron_name:\n",
    "        if neuron_to_cluster[i] != cluster:\n",
    "            continue\n",
    "        center = int(np.round(cortex_neuron[i]['soma_depth']*20))\n",
    "        distribution += func_gaussian(abs(distribution_index-center),size=kernel_size,sigma=kernel_sigma)\n",
    "    distribution /= np.sum(distribution)\n",
    "    plt.plot(np.convolve(distribution,sigma_kernel)[:1200],color=cluster_color_hex[cluster],label='cluster_'+str(cluster+1))\n",
    "plt.legend()\n",
    "plt.xlabel('Depth (\\u00B5m)',fontsize=15)\n",
    "plt.ylabel('Probability',fontsize=15,rotation=-90).set_verticalalignment('top')\n",
    "plt.title('Soma depth distribution',fontsize=17)\n",
    "plt.savefig(cluster_path+'cluster_depth_distribution.svg',format='svg')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analysis (only AId)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tmp = []\n",
    "for i in cortex_neuron:\n",
    "    if cortex_neuron[i]['cortex_subregion']=='d':\n",
    "        tmp.append(i)\n",
    "insular_neuron_name = tmp\n",
    "cluster_path = data_path+'cluster_allen/'\n",
    "create_dir(cluster_path)\n",
    "if_combined = False\n",
    "if if_combined:\n",
    "    neuron_axonal_length_by_region_log = []\n",
    "    for i in insular_neuron_name:\n",
    "        neuron_axonal_length_by_region_log.append([cortex_neuron[i]['axonal_projection_length_combined']])\n",
    "    neuron_axonal_length_by_region_log = np.array(neuron_axonal_length_by_region_log).reshape(-1,len(target_region))\n",
    "    for i in target_region_hierarchy['FT']:\n",
    "        neuron_axonal_length_by_region_log[:,target_region.index(i)] *= 2\n",
    "else:\n",
    "    neuron_axonal_length_by_region_log = []\n",
    "    for i in insular_neuron_name:\n",
    "        if cortex_neuron[i]['soma_hemisphere']=='left':\n",
    "            neuron_axonal_length_by_region_log.append([cortex_neuron[i]['axonal_projection_length_left'],cortex_neuron[i]['axonal_projection_length_right']])\n",
    "        else:\n",
    "            neuron_axonal_length_by_region_log.append([cortex_neuron[i]['axonal_projection_length_right'],cortex_neuron[i]['axonal_projection_length_left']])\n",
    "    neuron_axonal_length_by_region_log = np.array(neuron_axonal_length_by_region_log).reshape(-1,2*len(target_region))\n",
    "    for i in target_region_hierarchy['FT']:\n",
    "        neuron_axonal_length_by_region_log[:,target_region.index(i)] *= 2\n",
    "        neuron_axonal_length_by_region_log[:,target_region.index(i)+len(target_region)] *= 2\n",
    "    neuron_axonal_length_by_region_log = np.log(neuron_axonal_length_by_region_log+1)\n",
    "if not os.path.exists(cluster_path+'cluster_color_20.npy'):\n",
    "    cluster_color = np.array(distinctipy.get_colors(20,background_colors.tolist()))\n",
    "    np.save(cluster_path+'cluster_color_20.npy',cluster_color)\n",
    "else:\n",
    "    cluster_color = np.load(cluster_path+'cluster_color_20.npy')\n",
    "cluster_color_hex = []\n",
    "for i in cluster_color:\n",
    "    cluster_color_hex.append(rgb_to_hex(np.round(i*255).astype(int)))\n",
    "print('The colors of clustering:')\n",
    "distinctipy.color_swatch(cluster_color)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Clustering"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hierarchy.set_link_color_palette(cluster_color_hex)\n",
    "color_threhold = 200\n",
    "fig = plt.figure(figsize=(10,5))\n",
    "if not os.path.exists(cluster_path+'cluster_object2.npy'):\n",
    "    cluster_object = hierarchy.dendrogram(hierarchy.linkage(neuron_axonal_length_by_region_log, 'ward'),color_threshold=color_threhold,distance_sort=True)\n",
    "    np.save(cluster_path+'cluster_object.npy',cluster_object)\n",
    "else:\n",
    "    cluster_object = np.load(cluster_path+'cluster_object.npy',allow_pickle=True).item()\n",
    "plt.axhline(color_threhold, linestyle='dashed', color='#AAAAAA')\n",
    "plt.title(\"Hierachy of the clustering\",fontsize=20)\n",
    "plt.locator_params(axis='y', nbins=3)\n",
    "plt.xticks([],[])\n",
    "plt.xlabel('neurons = '+str(len(insular_neuron_name)),fontsize=15)\n",
    "plt.ylabel('Cluster distance',fontsize=15)\n",
    "#plt.savefig(cluster_path+'cluster_hierachy.svg',format='svg')\n",
    "plt.show()\n",
    "cluster_object_parse = [[],[]]\n",
    "neuron_to_cluster = {}\n",
    "for i,j in zip(cluster_object['leaves'],cluster_object['leaves_color_list']):\n",
    "    neuron = insular_neuron_name[i]\n",
    "    projection_matrix = neuron_axonal_length_by_region_log[i]\n",
    "    cluster_object_parse[0].append(cluster_color_hex.index(j))\n",
    "    cluster_object_parse[1].append(projection_matrix)\n",
    "    neuron_to_cluster[neuron] = cluster_color_hex.index(j)\n",
    "cluster_object_parse[1] = np.array(cluster_object_parse[1])\n",
    "cluster_object_parse_visualization = np.zeros(shape=(len(insular_neuron_name),1+len(target_region),3))\n",
    "cluster_object_parse_visualization[:,0,:] = np.array([cluster_color[i] for i in cluster_object_parse[0]])\n",
    "tmp = np.array([cluster_object_parse[1][:,:len(target_region)],cluster_object_parse[1][:,:len(target_region)]*0,cluster_object_parse[1][:,len(target_region):]]).transpose((1,2,0))\n",
    "tmp /= np.max(tmp)\n",
    "cluster_object_parse_visualization[:,1:,:] = tmp*1\n",
    "fig = plt.figure(figsize=(20,60))\n",
    "ax = fig.add_subplot()\n",
    "ax.matshow(cluster_object_parse_visualization.transpose((1,0,2)))\n",
    "_ = ax.set_xticks([],[])\n",
    "ax.set_xlabel('neurons = '+str(len(insular_neuron_name)),fontsize=fontsize_label)\n",
    "cluster_object_parse_visualization_yticks = [[0],['class']]\n",
    "cluster_object_parse_visualization_yticks[0] += [1+i for i in range(len(target_region))]\n",
    "cluster_object_parse_visualization_yticks[1] += target_region\n",
    "_ = ax.set_yticks(cluster_object_parse_visualization_yticks[0],cluster_object_parse_visualization_yticks[1],fontsize=fontsize_ticks)\n",
    "#plt.savefig(cluster_path+'cluster_projection_matrix.svg',format='svg')\n",
    "plt.show()\n",
    "tsne = TSNE(n_components=2, perplexity=35, random_state=0,n_iter=10000)\n",
    "X_tsne = tsne.fit_transform(neuron_axonal_length_by_region_log)\n",
    "tmp = {}\n",
    "for i,j in enumerate(insular_neuron_name):\n",
    "    _cluster = cluster_color_hex.index(cluster_object['leaves_color_list'][cluster_object['leaves'].index(i)])\n",
    "    if _cluster not in tmp:\n",
    "        tmp[_cluster] = []\n",
    "    tmp[_cluster].append([j,cortex_neuron[j]['soma_flatten_position'],X_tsne[i]])\n",
    "_flatmap = cv.imread('data/flatmapedge2.tif',1)\n",
    "fig = plt.figure(figsize=(10,10))\n",
    "flatmap = _flatmap*1\n",
    "for i in range(len(list(tmp.keys()))):\n",
    "    tmp2 = []\n",
    "    tmp3 = []\n",
    "    for j in tmp[i]:\n",
    "        if np.sum(abs(np.array(j[1])))!=0:\n",
    "            cv.circle(flatmap,[int(j[1][0])*2, int(j[1][1])*2],3,cluster_color[i]*255,-1)\n",
    "        tmp2.append(j[2])\n",
    "        tmp3.append([j[0][:-7],j[0][-7:]])\n",
    "    tmp2 = np.array(tmp2)\n",
    "    plt.scatter(tmp2[:,0], tmp2[:,1],c=cluster_color_hex[i],s=30,marker='o',label=str(i+1))\n",
    "    save_scene(tmp3,cluster_color[i],cluster_path+'cluster_'+str(i+1))\n",
    "plt.legend()\n",
    "plt.title(\"t-SNE Visualization of the clusters\",fontsize=fontsize_title)\n",
    "plt.axis('off')\n",
    "plt.savefig(cluster_path+'cluster_t-SNE.svg',format='svg')\n",
    "plt.show()\n",
    "fig = plt.figure(figsize=(10,10))\n",
    "plt.imshow(flatmap)\n",
    "plt.xlim(300,2700)\n",
    "plt.ylim(2700,300)\n",
    "plt.title(\"The soma distribution of the clusters\",fontsize=fontsize_title)\n",
    "plt.axis('off')\n",
    "#plt.savefig(cluster_path+'cluster_distribution.svg',format='svg')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "number_clusters = []\n",
    "for i in neuron_to_cluster:\n",
    "    if neuron_to_cluster[i] not in number_clusters:\n",
    "        number_clusters.append(neuron_to_cluster[i])\n",
    "number_clusters = len(number_clusters)\n",
    "sigma = 100\n",
    "sigma_kernel = [1/sigma for i in range(sigma)]\n",
    "for cluster in range(number_clusters):\n",
    "    distribution = np.zeros((1800,))\n",
    "    for i in insular_neuron_name:\n",
    "        if neuron_to_cluster[i] != cluster:\n",
    "            continue\n",
    "        distribution[int(np.round(cortex_neuron[i]['soma_depth']*20))] += 1\n",
    "    #distribution /= np.sum(distribution)\n",
    "    plt.plot(np.convolve(distribution,sigma_kernel)[:1200],color=cluster_color_hex[cluster],label=str(cluster+1))\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "insular_visualization_data_path = data_path+'visualization_allen/'\n",
    "create_dir(insular_visualization_data_path)\n",
    "visualization_sampling = 0.05\n",
    "neuronvis = nv.neuronVis()\n",
    "neuronvis.render.setBackgroundColor((1.0,1.0,1.0,1.0))\n",
    "control_console_output(False)\n",
    "for i in insular_neuron_name:\n",
    "    if np.random.uniform() < (1-visualization_sampling):\n",
    "        continue\n",
    "    neuron_tree = iondata.getNeuronTreeByID(i[:-7],i[-7:])\n",
    "    neuron_tree.somaRadius = 30\n",
    "    neuronvis.addNeuronTree(neuron_tree,somaColor=[0.059,0.616,0.345],dendriteColor=[0.258,0.52,0.957],axonColor=cluster_color[neuron_to_cluster[i]],somaHide=True,dendriteHide=True,mirrorToRight=True)\n",
    "    for j in ['anterior','dorsal','right']:\n",
    "        set_view(neuronvis,j)\n",
    "        neuronvis.render.savepng(insular_visualization_data_path+'all_'+j+'.png')\n",
    "control_console_output(True)\n",
    "neuronvis.render.closeWindow()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualization (each cluster)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "visualization_sampling = 0.05\n",
    "number_clusters = []\n",
    "for i in neuron_to_cluster:\n",
    "    if neuron_to_cluster[i] not in number_clusters:\n",
    "        number_clusters.append(neuron_to_cluster[i])\n",
    "number_clusters = len(number_clusters)\n",
    "neuronvis = nv.neuronVis()\n",
    "neuronvis.render.setBackgroundColor((1.0,1.0,1.0,1.0))\n",
    "control_console_output(False)\n",
    "for cluster in range(number_clusters):\n",
    "    neuronvis.clear()\n",
    "    for i in insular_neuron_name:\n",
    "        if neuron_to_cluster[i] != cluster:\n",
    "            continue\n",
    "        if np.random.uniform() < (1-visualization_sampling):\n",
    "            continue\n",
    "        neuron_tree = iondata.getNeuronTreeByID(i[:-7],i[-7:])\n",
    "        neuron_tree.somaRadius = 30\n",
    "        neuronvis.addNeuronTree(neuron_tree,somaColor=[0.059,0.616,0.345],dendriteColor=[0.258,0.52,0.957],axonColor=cluster_color[neuron_to_cluster[i]],somaHide=True,dendriteHide=True,mirrorToRight=True)\n",
    "    for j in ['anterior','dorsal','right']:\n",
    "        set_view(neuronvis,j)\n",
    "        neuronvis.render.savepng(insular_visualization_data_path+'cluster_'+str(cluster+1)+'_'+j+'.png')\n",
    "control_console_output(True)\n",
    "neuronvis.render.closeWindow()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## GPe analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "GPe_projection_matrix = []\n",
    "for i in insular_neuron_name:\n",
    "        GPe_projection_matrix.append(get_axonal_statics_by_regions(i,['GPe'],'projectleftregion')+get_axonal_statics_by_regions(i,['GPe'],'projectrightregion'))\n",
    "GPe_projection_matrix = np.array(GPe_projection_matrix)\n",
    "GPe_projection_matrix_log = np.log(GPe_projection_matrix+1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tmp = GPe_projection_matrix_log/np.max(GPe_projection_matrix_log)\n",
    "tmp2 = []\n",
    "for i in cluster_object['leaves']:\n",
    "    j = tmp[i]\n",
    "    tmp2.append([j[0],0,j[1]])\n",
    "tmp2 = np.array(tmp2).reshape(1,-1,3)\n",
    "plt.matshow(tmp2)\n",
    "plt.axis('off')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "k = np.power(GPe_projection_matrix,0.25)\n",
    "tmp = k/np.max(k)\n",
    "tmp2 = []\n",
    "for i in cluster_object['leaves']:\n",
    "    j = tmp[i]\n",
    "    tmp2.append([j[0],0,j[1]])\n",
    "tmp2 = np.array(tmp2).reshape(1,-1,3)\n",
    "plt.matshow(tmp2)\n",
    "plt.axis('off')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "GPe_projection_matrix = []\n",
    "for i in insular_neuron_name:\n",
    "        GPe_projection_matrix.append(get_axonal_statics_by_regions(i,['GPe'],'terminalleftregion')+get_axonal_statics_by_regions(i,['GPe'],'terminalrightregion'))\n",
    "GPe_projection_matrix = np.array(GPe_projection_matrix)\n",
    "GPe_projection_matrix_log = np.log(GPe_projection_matrix+1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "k = np.power(GPe_projection_matrix,0.5)\n",
    "tmp = k/np.max(k)\n",
    "tmp2 = []\n",
    "for i in cluster_object['leaves']:\n",
    "    j = tmp[i]\n",
    "    tmp2.append([j[0],0,j[1]])\n",
    "tmp2 = np.array(tmp2).reshape(1,-1,3)\n",
    "plt.matshow(tmp2)\n",
    "plt.axis('off')\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tmp = {}\n",
    "for i,j in zip(insular_neuron_name,GPe_projection_matrix):\n",
    "    if np.sum(j)==0:\n",
    "        continue\n",
    "    cluster = neuron_to_cluster[i]\n",
    "    if cluster not in tmp:\n",
    "        tmp[cluster] = []\n",
    "    tmp[cluster].append([i[:-7],i[-7:]])\n",
    "for i in tmp:\n",
    "    save_scene(tmp[i],cluster_color[i],cluster_path+'GPe_cluster_'+str(i+1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Xiao version analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#visualization of each neuron by xiao's segmentation\n",
    "insular_visualization_data_path = data_path+'visualization/'\n",
    "create_dir(insular_visualization_data_path)\n",
    "insular_visualization_path = {}\n",
    "for i in region_transform_region:\n",
    "    insular_visualization_path[i] = insular_visualization_data_path+i+'/'\n",
    "    create_dir(insular_visualization_path[i])\n",
    "neuronvis = nv.neuronVis()\n",
    "neuronvis.render.setBackgroundColor((1.0,1.0,1.0,1.0))\n",
    "for i in insular_neuron_name:\n",
    "    neuronvis.clear()\n",
    "    neuron_tree = iondata.getNeuronTreeByID(i[:-7],i[-7:])\n",
    "    neuron_tree.somaRadius = 30\n",
    "    neuronvis.addNeuronTree(neuron_tree,somaColor=[0.059,0.616,0.345],dendriteColor=[0.258,0.52,0.957],axonColor=[0.859,0.267,0.216],mirrorToRight=True)\n",
    "    bregma = round(-(cortex_neuron[i]['soma_position'][0]-5472)/1000,2)\n",
    "    for j in ['dorsal','right']:\n",
    "        set_view(neuronvis,j)\n",
    "        neuronvis.render.savepng(insular_visualization_path[region_transform[i]]+i[:-4]+'_'+str(bregma)+'_'+j+'.png')\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cluster_path = data_path+'cluster_xiao/'\n",
    "create_dir(cluster_path)\n",
    "if_combined = False\n",
    "if if_combined:\n",
    "    neuron_axonal_length_by_region_log = []\n",
    "    for i in insular_neuron_name:\n",
    "        neuron_axonal_length_by_region_log.append([cortex_neuron[i]['axonal_projection_length_combined']])\n",
    "    neuron_axonal_length_by_region_log = np.array(neuron_axonal_length_by_region_log).reshape(-1,len(target_region))\n",
    "    for i in target_region_hierarchy['FT']:\n",
    "        neuron_axonal_length_by_region_log[:,target_region.index(i)] *= 2\n",
    "else:\n",
    "    neuron_axonal_length_by_region_log = []\n",
    "    for i in insular_neuron_name:\n",
    "        if cortex_neuron[i]['soma_hemisphere']=='left':\n",
    "            neuron_axonal_length_by_region_log.append([cortex_neuron[i]['axonal_projection_length_left'],cortex_neuron[i]['axonal_projection_length_right']])\n",
    "        else:\n",
    "            neuron_axonal_length_by_region_log.append([cortex_neuron[i]['axonal_projection_length_right'],cortex_neuron[i]['axonal_projection_length_left']])\n",
    "    neuron_axonal_length_by_region_log = np.array(neuron_axonal_length_by_region_log).reshape(-1,2*len(target_region))\n",
    "    for i in target_region_hierarchy['FT']:\n",
    "        neuron_axonal_length_by_region_log[:,target_region.index(i)] *= 2\n",
    "        neuron_axonal_length_by_region_log[:,target_region.index(i)+len(target_region)] *= 2\n",
    "    neuron_axonal_length_by_region_log = np.log(neuron_axonal_length_by_region_log+1)\n",
    "if not os.path.exists(cluster_path+'cluster_color_20.npy'):\n",
    "    cluster_color = np.array(distinctipy.get_colors(20,background_colors.tolist()))\n",
    "    np.save(cluster_path+'cluster_color_20.npy',cluster_color)\n",
    "else:\n",
    "    cluster_color = np.load(cluster_path+'cluster_color_20.npy')\n",
    "cluster_color_hex = []\n",
    "for i in cluster_color:\n",
    "    cluster_color_hex.append(rgb_to_hex(np.round(i*255).astype(int)))\n",
    "print('The colors of clustering:')\n",
    "distinctipy.color_swatch(cluster_color_hex)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#overall information of the clustering\n",
    "hierarchy.set_link_color_palette(cluster_color_hex)\n",
    "color_threhold = 200\n",
    "fig = plt.figure(figsize=(10,5))\n",
    "if not os.path.exists(cluster_path+'cluster_object2.npy'):\n",
    "    cluster_object = hierarchy.dendrogram(hierarchy.linkage(neuron_axonal_length_by_region_log, 'ward'),color_threshold=color_threhold,distance_sort=True)\n",
    "    np.save(cluster_path+'cluster_object.npy',cluster_object)\n",
    "else:\n",
    "    cluster_object = np.load(cluster_path+'cluster_object.npy',allow_pickle=True).item()\n",
    "plt.axhline(color_threhold, linestyle='dashed', color='#AAAAAA')\n",
    "plt.title(\"Hierachy of the clustering\",fontsize=20)\n",
    "plt.locator_params(axis='y', nbins=3)\n",
    "plt.xticks([],[])\n",
    "plt.xlabel('neurons = '+str(len(insular_neuron_name)),fontsize=15)\n",
    "plt.ylabel('Cluster distance',fontsize=15)\n",
    "#plt.savefig(cluster_path+'cluster_hierachy.svg',format='svg')\n",
    "plt.show()\n",
    "cluster_object_parse = [[],[]]\n",
    "for i,j in zip(cluster_object['leaves'],cluster_object['leaves_color_list']):\n",
    "    neuron = insular_neuron_name[i]\n",
    "    projection_matrix = neuron_axonal_length_by_region_log[i]\n",
    "    cluster_object_parse[0].append(cluster_color_hex.index(j))\n",
    "    cluster_object_parse[1].append(projection_matrix)\n",
    "cluster_object_parse[1] = np.array(cluster_object_parse[1])\n",
    "cluster_object_parse_visualization = np.zeros(shape=(len(insular_neuron_name),1+len(target_region),3))\n",
    "cluster_object_parse_visualization[:,0,:] = np.array([cluster_color[i] for i in cluster_object_parse[0]])\n",
    "tmp = np.array([cluster_object_parse[1][:,:len(target_region)],cluster_object_parse[1][:,:len(target_region)]*0,cluster_object_parse[1][:,len(target_region):]]).transpose((1,2,0))\n",
    "tmp /= np.max(tmp)\n",
    "cluster_object_parse_visualization[:,1:,:] = tmp*1\n",
    "fig = plt.figure(figsize=(20,60))\n",
    "ax = fig.add_subplot()\n",
    "ax.matshow(cluster_object_parse_visualization.transpose((1,0,2)))\n",
    "_ = ax.set_xticks([],[])\n",
    "ax.set_xlabel('neurons = '+str(len(insular_neuron_name)),fontsize=fontsize_label)\n",
    "cluster_object_parse_visualization_yticks = [[0],['class']]\n",
    "cluster_object_parse_visualization_yticks[0] += [1+i for i in range(len(target_region))]\n",
    "cluster_object_parse_visualization_yticks[1] += target_region\n",
    "_ = ax.set_yticks(cluster_object_parse_visualization_yticks[0],cluster_object_parse_visualization_yticks[1],fontsize=fontsize_ticks)\n",
    "#plt.savefig(cluster_path+'cluster_projection_matrix.svg',format='svg')\n",
    "plt.show()\n",
    "tsne = TSNE(n_components=2, perplexity=35, random_state=0,n_iter=10000)\n",
    "X_tsne = tsne.fit_transform(neuron_axonal_length_by_region_log)\n",
    "tmp = {}\n",
    "for i,j in enumerate(insular_neuron_name):\n",
    "    _cluster = cluster_color_hex.index(cluster_object['leaves_color_list'][cluster_object['leaves'].index(i)])\n",
    "    if _cluster not in tmp:\n",
    "        tmp[_cluster] = []\n",
    "    tmp[_cluster].append([j,cortex_neuron[j]['soma_flatten_position'],X_tsne[i]])\n",
    "_flatmap = cv.imread('data/flatmapedge2.tif',1)\n",
    "fig = plt.figure(figsize=(10,10))\n",
    "flatmap = _flatmap*1\n",
    "for i in range(len(list(tmp.keys()))):\n",
    "    tmp2 = []\n",
    "    tmp3 = []\n",
    "    for j in tmp[i]:\n",
    "        if np.sum(abs(np.array(j[1])))!=0:\n",
    "            cv.circle(flatmap,[int(j[1][0])*2, int(j[1][1])*2],7,cluster_color[i]*255,-1)\n",
    "        tmp2.append(j[2])\n",
    "        tmp3.append([j[0][:-7],j[0][-7:]])\n",
    "    tmp2 = np.array(tmp2)\n",
    "    plt.scatter(tmp2[:,0], tmp2[:,1],c=cluster_color_hex[i],s=30,marker='o',label=str(i+1))\n",
    "    save_scene(tmp3,cluster_color[i],cluster_path+'cluster_'+str(i+1))\n",
    "plt.legend()\n",
    "plt.title(\"t-SNE Visualization of the clusters\",fontsize=fontsize_title)\n",
    "plt.axis('off')\n",
    "plt.savefig(cluster_path+'cluster_t-SNE.svg',format='svg')\n",
    "plt.show()\n",
    "fig = plt.figure(figsize=(10,10))\n",
    "plt.imshow(flatmap)\n",
    "plt.xlim(300,2700)\n",
    "plt.ylim(2700,300)\n",
    "plt.title(\"The soma distribution of the clusters\",fontsize=fontsize_title)\n",
    "plt.axis('off')\n",
    "#plt.savefig(cluster_path+'cluster_distribution.svg',format='svg')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mode = 1\n",
    "angle_threshold = 45\n",
    "length_ratio_threshold = 0.8\n",
    "length_threshold = -1\n",
    "return_nonshaft = False\n",
    "mirror = True\n",
    "neuronvis = nv.neuronVis()\n",
    "neuronvis.render.setBackgroundColor((1.0,1.0,1.0,1.0))\n",
    "neuronvis_shaft = nv.neuronVis()\n",
    "neuronvis_shaft.render.setBackgroundColor((1.0,1.0,1.0,1.0))\n",
    "for i in range(len(list(tmp.keys()))):\n",
    "    final_lines = []\n",
    "    final_colors = []\n",
    "    final_lines2 = []\n",
    "    final_colors2 = []\n",
    "    neuronvis.clear()\n",
    "    neuronvis_shaft.clear()\n",
    "    tmp_distribution = [0 for i in region_transform_region]\n",
    "    flatmap = _flatmap*1\n",
    "    control_console_output(False)\n",
    "    for j in tmp[i]:\n",
    "        if np.sum(abs(np.array(j[1])))!=0:\n",
    "            cv.circle(flatmap,[int(j[1][0])*2, int(j[1][1])*2],15,cluster_color[i]*255,-1)\n",
    "        tmp_distribution[region_transform_region.index(region_transform[j[0]])]  += 1#need to change\n",
    "        neuron_tree = iondata.getNeuronTreeByID(j[0][:-7],j[0][-7:])\n",
    "        neuronvis.addNeuronTree(neuron_tree,axonColor=cluster_color[i],somaHide=True,dendriteHide=True,mirrorToRight=True,isLine=True)\n",
    "        [[lines,lines_color],[lines2,lines_color2]] = get_shaft(neuron_tree,mode=mode,return_nonshaft=return_nonshaft,angle_threshold=45,length_ratio_threshold=length_ratio_threshold,length_threshold=length_threshold,mirror=mirror)\n",
    "        final_lines += lines\n",
    "        final_colors += lines_color\n",
    "        final_lines2 += lines2\n",
    "        final_colors2 += lines_color2\n",
    "    control_console_output(True)\n",
    "    fig = plt.figure(figsize=(25,10))\n",
    "    grid = fig.add_gridspec(1, 2, width_ratios=[1, 1.3])\n",
    "    ax1 = fig.add_subplot(grid[0])\n",
    "    ax1.imshow(flatmap)\n",
    "    ax1.set_xlim(300,2700)\n",
    "    ax1.set_ylim(2700,300)\n",
    "    ax1.set_axis_off()\n",
    "    ax2 = fig.add_subplot(grid[1])\n",
    "    ax2.bar(range(len(tmp_distribution)),tmp_distribution)\n",
    "    ax2.set_xticks(range(len(tmp_distribution)),region_transform_region,rotation=-90)\n",
    "    ax2.set_ylabel('Num. of neurons',fontsize=15)\n",
    "    plt.title(\"The soma distribution of the clusters \"+str(i+1),fontsize=20)\n",
    "    plt.savefig(cluster_path+'cluster_distribution_'+str(i+1)+'.svg',format='svg')\n",
    "    plt.show()\n",
    "    draw_lines(neuronvis_shaft,final_lines,final_colors)\n",
    "    if return_nonshaft:\n",
    "        draw_lines(neuronvis_shaft,final_lines2,final_colors2)\n",
    "    for j in ['dorsal','right']:\n",
    "        set_view(neuronvis,j)\n",
    "        neuronvis.render.savepng(cluster_path+'cluster_'+str(i+1)+'_axon_'+j+'.png')\n",
    "    for j in ['dorsal','right']:\n",
    "        set_view(neuronvis_shaft,j)\n",
    "        neuronvis_shaft.render.savepng(cluster_path+'cluster_'+str(i+1)+'_axon_shaft_'+j+'.png')\n",
    "neuronvis.render.closeWindow()\n",
    "neuronvis_shaft.render.closeWindow()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def hex_to_rgb(hex_code):\n",
    "    # Remove the '#' symbol if present\n",
    "    hex_code = hex_code.lstrip('#')\n",
    "    \n",
    "    # Convert the hexadecimal string to integer values\n",
    "    red = int(hex_code[0:2], 16)\n",
    "    green = int(hex_code[2:4], 16)\n",
    "    blue = int(hex_code[4:6], 16)\n",
    "    \n",
    "    return red, green, blue\n",
    "\n",
    "depth = []\n",
    "position = []\n",
    "for i in insular_neuron_name:\n",
    "    depth.append([cortex_neuron[i]['soma_depth']])\n",
    "    position.append(cortex_neuron[i]['soma_flatten_position'])\n",
    "position_color = position_color_mapping(np.array(position),2)*1.5\n",
    "depth_color = position_color_mapping(np.array(depth),1)*2.5\n",
    "_flatmap = cv.imread('data/flatmapedge2.tif',1)\n",
    "fig = plt.figure(figsize=(10,10))\n",
    "flatmap = _flatmap*0+255\n",
    "_max_x = np.max(abs(X_tsne[:,0]))+10\n",
    "_max_y = np.max(abs(X_tsne[:,1]))+10\n",
    "for i,j in zip(X_tsne,insular_neuron_name):\n",
    "    c = region_transform_color[region_transform[j]]\n",
    "    if np.sum(abs(np.array(i)))!=0:\n",
    "        cv.circle(flatmap,[int(i[0]/_max_x*1500+1500),int(i[1]/_max_y*1500+1500)],15,hex_to_rgb(c),-1)\n",
    "plt.matshow(flatmap[::-1,...])\n",
    "plt.xlim(300,2700)\n",
    "plt.ylim(2700,300)\n",
    "plt.title(\"t-SNE Visualization (corlored by a-m-p-v-g)\",fontsize=fontsize_title*0.4)\n",
    "plt.axis('off')\n",
    "plt.show()\n",
    "fig = plt.figure(figsize=(10,10))\n",
    "flatmap = _flatmap*1\n",
    "for i,j in zip(position,position_color):\n",
    "    if np.sum(abs(np.array(i)))!=0:\n",
    "        cv.circle(flatmap,[int(i[0])*2, int(i[1])*2],7,j*255,-1)\n",
    "plt.matshow(flatmap)\n",
    "plt.xlim(300,2700)\n",
    "plt.ylim(2700,300)\n",
    "plt.title(\"t-SNE Visualization of the clusters\",fontsize=fontsize_title*0.4)\n",
    "plt.axis('off')\n",
    "plt.show()\n",
    "fig = plt.figure(figsize=(10,10))\n",
    "flatmap = _flatmap*0+255\n",
    "_max_x = np.max(abs(X_tsne[:,0]))+10\n",
    "_max_y = np.max(abs(X_tsne[:,1]))+10\n",
    "for i,j in zip(X_tsne,position_color):\n",
    "    if np.sum(abs(np.array(i)))!=0:\n",
    "        cv.circle(flatmap,[int(i[0]/_max_x*1500+1500),int(i[1]/_max_y*1500+1500)],15,j*255,-1)\n",
    "plt.matshow(flatmap[::-1,...])\n",
    "plt.xlim(300,2700)\n",
    "plt.ylim(2700,300)\n",
    "plt.title(\"t-SNE Visualization (corlored by position)\",fontsize=fontsize_title*0.4)\n",
    "plt.axis('off')\n",
    "plt.show()\n",
    "fig = plt.figure(figsize=(10,10))\n",
    "flatmap = _flatmap*0+255\n",
    "_max_x = np.max(abs(X_tsne[:,0]))+10\n",
    "_max_y = np.max(abs(X_tsne[:,1]))+10\n",
    "for i,j in zip(X_tsne,depth_color):\n",
    "    if np.sum(abs(np.array(i)))!=0:\n",
    "        cv.circle(flatmap,[int(i[0]/_max_x*1500+1500),int(i[1]/_max_y*1500+1500)],15,j*255,-1)\n",
    "plt.matshow(flatmap[::-1,...])\n",
    "plt.xlim(300,2700)\n",
    "plt.ylim(2700,300)\n",
    "plt.title(\"t-SNE Visualization (corlored by soma depth)\",fontsize=fontsize_title*0.4)\n",
    "plt.axis('off')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tmp = seaborn.clustermap(np.log(projection_matrix+1).T,row_cluster=False)\n",
    "#tmp.ax_heatmap.set_yticklabels(target_region)\n",
    "plt.show()\n",
    "projection_matrix_visualization = []\n",
    "projection_matrix_indices_visualization = []\n",
    "row_indices = tmp.data2d.index.tolist()\n",
    "column_indices = tmp.data2d.columns.tolist()\n",
    "for i in column_indices:\n",
    "    projection_matrix_indices_visualization.append(index_bar[i])\n",
    "projection_matrix_indices_visualization = np.array(projection_matrix_indices_visualization)\n",
    "_max = np.max(projection_matrix)\n",
    "for i in column_indices:\n",
    "    projection_matrix_visualization.append(projection_matrix[i])\n",
    "projection_matrix_visualization = np.array(projection_matrix_visualization)\n",
    "plt.matshow(projection_matrix_indices_visualization.reshape(1,-1,3)*255)\n",
    "plt.axis('off')\n",
    "plt.show()\n",
    "icpt = []\n",
    "icpt_color = {'IT':[1,0,0],'CT':[0,1,0],'PT':[0,0,1],'Others':[1,1,0]}\n",
    "for i in column_indices:\n",
    "    icpt.append(icpt_color[cortex_neuron[index_name[i]]['type']])\n",
    "icpt = np.array(icpt)\n",
    "plt.matshow(icpt.reshape(1,-1,3)*255)\n",
    "plt.axis('off')\n",
    "plt.show()\n",
    "tmp_color = np.array(distinctipy.get_colors(len(target_region_hierarchy.keys()),background_colors.tolist()))\n",
    "tmp = []\n",
    "for ii,i in enumerate(target_region_hierarchy):\n",
    "    for j in target_region_hierarchy[i]:\n",
    "        tmp.append(tmp_color[ii])\n",
    "tmp_r = np.array(tmp)\n",
    "plt.matshow(tmp_r.reshape(1,-1,3))\n",
    "plt.axis('off')\n",
    "plt.show()\n",
    "fig = plt.figure(figsize=(15,15))\n",
    "ax = fig.add_subplot(1,1,1)\n",
    "tmp = ax.matshow(np.log(projection_matrix_visualization+1).T,cmap='hot')\n",
    "ax.set_xticks([],[])\n",
    "ax.set_xlabel('Num. of neuron =  '+str(projection_matrix_visualization.shape[0]),fontsize=fontsize_label)\n",
    "ax.set_yticks(range(len(target_region)),target_region,fontsize=fontsize_ticks/3)\n",
    "ax.set_ylabel('Brain region',fontsize=fontsize_label)\n",
    "plt.colorbar(tmp)\n",
    "plt.show()\n",
    "icpt_region = {'AIa':['AId','AIv'],'AIp':['AIp'],'GI':['VISC','GU']}\n",
    "icpt = {'AIa':[0,0,0,0],'AIp':[0,0,0,0],'GI':[0,0,0,0]}\n",
    "icpt_index = ['IT','CT','PT','Others']\n",
    "icpt_color = {'IT':[1,0,0],'CT':[0,1,0],'PT':[0,0,1],'Others':[1,1,0]}\n",
    "icpt_color_hex = []\n",
    "for i in icpt_color:\n",
    "    icpt_color_hex.append(rgb_to_hex(np.array(icpt_color[i])*255))\n",
    "for i in cortex_neuron:\n",
    "    tmp = cortex_neuron[i]['cortex_region']+cortex_neuron[i]['cortex_subregion']\n",
    "    for j in icpt_region:\n",
    "        if tmp in icpt_region[j]:\n",
    "            icpt[j][icpt_index.index(cortex_neuron[i]['type'])] += 1\n",
    "fig = plt.figure(figsize=(15,5))\n",
    "for ii,i in enumerate(icpt):\n",
    "    ax = fig.add_subplot(1,3,ii+1)\n",
    "    ax.pie(icpt[i][:3],labels=icpt_index[:3],colors=icpt_color_hex[:3],autopct='%1.1f%%')\n",
    "    ax.set_title(i,fontsize=fontsize_title)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#dendrites\n",
    "for i in cortex_neuron:\n",
    "    if 1:\n",
    "        j = i\n",
    "        if np.random.rand()<0.8:\n",
    "            continue\n",
    "        neuron_tree = iondata.getNeuronTreeByID(j[:-7],j[-7:])\n",
    "        mark = 0\n",
    "        for edge in neuron_tree.edges:\n",
    "            if edge.data[0].type!=3:\n",
    "                continue\n",
    "            mark = 1\n",
    "            tmp = []\n",
    "            for point in edge.data[::20]:\n",
    "                tmp.append(get_mirrored_point(point.xyz))\n",
    "            tmp.append(get_mirrored_point(edge.data[-1].xyz))\n",
    "            tmp = compute_flatten_position_and_depth(np.array(tmp))\n",
    "            plt.plot([i[3][0] for i in tmp],[i[0] for i in tmp],'g',linewidth=0.5)\n",
    "        if mark==0:\n",
    "            continue\n",
    "        plt.scatter([cortex_neuron[i]['soma_flatten_position'][0]],[cortex_neuron[i]['soma_depth']],color='r',s=1)\n",
    "        plt.title(cortex_neuron[i]['cortex_region']+cortex_neuron[i]['cortex_subregion'])\n",
    "        plt.ylim(60,0)\n",
    "        plt.axis('off')\n",
    "        plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#figure 1\n",
    "tmp = {'AIa':['AId','AIv'],'AIp':['AIp'],'GI':['VISC','GU']}\n",
    "tmp_color = {'AIa':'#DB4537','AIp':'#4286F3','GI':'#55AF7B'}\n",
    "allen_to_xiao = {}\n",
    "for i in tmp:\n",
    "    for j in tmp[i]:\n",
    "        allen_to_xiao[j] = i\n",
    "tmp = {}\n",
    "for i in cortex_neuron:\n",
    "    region = allen_to_xiao[cortex_neuron[i]['cortex_region']+cortex_neuron[i]['cortex_subregion']]\n",
    "    if region not in tmp:\n",
    "        tmp[region] = {}\n",
    "    layer = cortex_neuron[i]['cortex_layer']\n",
    "    if layer not in tmp[region]:\n",
    "        tmp[region][layer] = {}\n",
    "    type = cortex_neuron[i]['type']\n",
    "    if type not in tmp[region][layer]:\n",
    "        tmp[region][layer][type] = 0.0\n",
    "    tmp[region][layer][type] += 1.0\n",
    "fig = plt.figure(figsize=(15,3))\n",
    "ax = fig.add_subplot(1,3,1)\n",
    "for region_index,region in enumerate(['AIa','AIp','GI']):\n",
    "    tmp2 = []\n",
    "    for layer in ['1','2/3','5','6']:\n",
    "        tmp2.append(0.0)\n",
    "        for type in ['IT','CT','PT']:\n",
    "            try:\n",
    "                tmp2[-1] += tmp[region][layer][type]\n",
    "            except:\n",
    "                pass\n",
    "    tmp2 = np.array(tmp2)\n",
    "    tmp2 /= np.sum(tmp2)\n",
    "    tmp2 *= 100\n",
    "    ax.bar([i*2+region_index/2 for i in range(4)],tmp2,width=0.45,label=region,color=tmp_color[region],capstyle='round')\n",
    "ax.set_xticks([i*2+0.5 for i in range(4)],['1','2/3','5','6'],fontsize=fontsize_ticks*0.6)\n",
    "ax.set_yticks([0,50,100],[0,50,100],fontsize=fontsize_ticks*0.6)\n",
    "ax.set_ylim([0,100])\n",
    "ax.set_xlabel('Layer',fontsize=fontsize_label*0.6)\n",
    "ax.set_ylabel('percentage (%)',fontsize=fontsize_label*0.6,rotation=90)\n",
    "ax.set_title('Soma distribution',fontsize=fontsize_title*0.6)\n",
    "ax.legend()\n",
    "ax = fig.add_subplot(1,3,3)\n",
    "for region_index,region in enumerate(['AIa','AIp','GI']):\n",
    "    tmp2 = []\n",
    "    for type in ['IT','CT','PT']:\n",
    "        tmp2.append(0.0)\n",
    "        for layer in ['1','2/3','5','6']:\n",
    "            try:\n",
    "                tmp2[-1] += tmp[region][layer][type]\n",
    "            except:\n",
    "                pass\n",
    "    tmp2 = np.array(tmp2)\n",
    "    tmp2 /= np.sum(tmp2)\n",
    "    tmp2 *= 100\n",
    "    ax.bar([i*2+region_index/2 for i in range(3)],tmp2,width=0.45,label=region,color=tmp_color[region],capstyle='round')\n",
    "ax.set_xticks([i*2+0.5 for i in range(3)],['IT','CT','PT'],fontsize=fontsize_ticks*0.6)\n",
    "ax.set_yticks([0,50,100],[0,50,100],fontsize=fontsize_ticks*0.6)\n",
    "ax.set_ylim([0,100])\n",
    "ax.set_xlabel('Type',fontsize=fontsize_label*0.6)\n",
    "ax.set_ylabel('percentage (%)',fontsize=fontsize_label*0.6,rotation=90)\n",
    "ax.set_title('Soma distribution',fontsize=fontsize_title*0.6)\n",
    "ax.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#figure 3\n",
    "tmp = {'AIa':['AId','AIv'],'AIp':['AIp'],'GI':['VISC','GU']}\n",
    "tmp_color = {'AIa':'#DB4537','AIp':'#4286F3','GI':'#55AF7B'}\n",
    "allen_to_xiao = {}\n",
    "for i in tmp:\n",
    "    for j in tmp[i]:\n",
    "        allen_to_xiao[j] = i\n",
    "tmp = [[i,j] for i, j in sorted(enumerate(np.mean(projection_matrix,axis=0)), key=lambda x: x[1], reverse=True)]\n",
    "_cortical_region_data = []\n",
    "_cortical_region_data_label = []\n",
    "_subcortical_region_data = []\n",
    "_subcortical_region_data_label = []\n",
    "_cortical_region = target_region_hierarchy['Cortex']+target_region_hierarchy['OLF']+target_region_hierarchy['HPF']+target_region_hierarchy['Cortex-sp']\n",
    "for i in range(20):\n",
    "    region = tmp[i][0]\n",
    "    if target_region[region] in _cortical_region:\n",
    "        tmp2 = _cortical_region_data\n",
    "        tmp3 = _cortical_region_data_label\n",
    "    else:\n",
    "        tmp2 = _subcortical_region_data\n",
    "        tmp3 = _subcortical_region_data_label\n",
    "    tmp2.append([0,0,0])\n",
    "    tmp3.append(target_region[region])\n",
    "    for j in cortex_neuron:\n",
    "        cortex_neuron_region_index = ['AIa','AIp','GI'].index(allen_to_xiao[cortex_neuron[j]['cortex_region']+cortex_neuron[j]['cortex_subregion']])\n",
    "        tmp2[-1][cortex_neuron_region_index] += cortex_neuron[j]['axonal_projection_length_combined'][region]\n",
    "fig = plt.figure(figsize=(20,4))\n",
    "ax = fig.add_subplot(1,1,1)\n",
    "for i,j in zip(range(3),['AIa','AIp','GI']):\n",
    "    #ax.bar([k*2+i/2 for k in range(20)],np.log(np.array(_cortical_region_data+_subcortical_region_data)[:,i]+1),width=0.45,label=j,color=tmp_color[j],capstyle='round')\n",
    "    ax.bar([k*2+i/2 for k in range(20)],np.power(np.array(_cortical_region_data+_subcortical_region_data)[:,i],0.5),width=0.45,label=j,color=tmp_color[j],capstyle='round')\n",
    "ax.set_xticks([i*2+0.5 for i in range(20)],_cortical_region_data_label+_subcortical_region_data_label,fontsize=fontsize_ticks*0.6)\n",
    "#ax.set_yticks([0,50,100],[0,50,100],fontsize=fontsize_ticks*0.6)\n",
    "ax.set_xlabel('Target region',fontsize=fontsize_label*0.6)\n",
    "ax.set_ylabel('Axonal length ($\\sqrt{\\u00B5m}$)',fontsize=fontsize_label*0.6,rotation=90)\n",
    "ax.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#figure 4\n",
    "iso_mask = get_region_mask('Isocortex')\n",
    "unit = 10\n",
    "sampling = 40\n",
    "projection_matrix_size = np.array(iso_mask.shape).astype(int)//unit\n",
    "slected_neuron_region = {'AIa':['AId','AIv'],'AIp':['AIp'],'GI':['GU','VISC']}\n",
    "slected_neuron_region_index = []\n",
    "projection_matrix = {}\n",
    "for i in slected_neuron_region:\n",
    "    slected_neuron_region_index.append(i)\n",
    "    projection_matrix[i] = {}\n",
    "    projection_matrix[i]['matrix'] = np.zeros(projection_matrix_size[:2])\n",
    "    projection_matrix[i]['num'] = 0\n",
    "for i in tool_bar(cortex_neuron):\n",
    "    index = ''\n",
    "    for j in slected_neuron_region:\n",
    "        if cortex_neuron[i]['cortex_region']+cortex_neuron[i]['cortex_subregion'] in slected_neuron_region[j]:\n",
    "            index = j\n",
    "            break\n",
    "    if index=='':\n",
    "        continue\n",
    "    control_console_output(0)\n",
    "    neuron_tree = iondata.getNeuronTreeByID(i[:-7],i[-7:])\n",
    "    for point in neuron_tree.points[::sampling]:\n",
    "        if point.type==3:\n",
    "            continue\n",
    "        tmp = get_mirrored_point(point.xyz)\n",
    "        if if_inside_mask(iso_mask,tmp/10):\n",
    "            projection_matrix[index]['matrix'][int(tmp[0]/10/unit),int(tmp[1]/10/unit)] += 1\n",
    "    projection_matrix[index]['num'] += 1\n",
    "    control_console_output(1)\n",
    "np.save('projection_matrix.npy',projection_matrix)\n",
    "tmp = []\n",
    "for i in projection_matrix:\n",
    "    tmp.append(projection_matrix[i]['matrix'].T/np.max(projection_matrix[i]['matrix']))\n",
    "c = tmp[2]\n",
    "tmp[2] = tmp[1]\n",
    "tmp[1] = c\n",
    "tmp = np.array(tmp).transpose(1,2,0)\n",
    "tmp = np.log(tmp+1)\n",
    "tmp /= np.max(tmp)\n",
    "fig = plt.figure(figsize=(10, 8))\n",
    "gs = fig.add_gridspec(2, 2,  width_ratios=(4.25, 0.75), height_ratios=(0.75, 3.25),left=0.1, right=0.9, bottom=0.1, top=0.9,wspace=0.05, hspace=0.05)\n",
    "ax = fig.add_subplot(gs[1, 0])\n",
    "ax_histx = fig.add_subplot(gs[0, 0], sharex=ax)\n",
    "ax_histy = fig.add_subplot(gs[1, 1], sharey=ax)\n",
    "#tmp = np.where(tmp==0,1,tmp)\n",
    "#tmp += np.where(np.sum(tmp,axis=-1)==0,1,0).reshape(tmp.shape[0],tmp.shape[1],1)\n",
    "ax.matshow(tmp[:,10:110])\n",
    "for i,j,k in zip(projection_matrix,['r','b','g'],['AIa','AIp','GI']):\n",
    "    tmp_distribution = projection_matrix[i]['matrix'][10:110]\n",
    "    tmp_distribution_x = np.sum(tmp_distribution,axis=1)\n",
    "    tmp_distribution_x /= np.linalg.norm(tmp_distribution_x)\n",
    "    tmp_distribution_y = np.sum(tmp_distribution,axis=0)\n",
    "    tmp_distribution_y /= np.linalg.norm(tmp_distribution_y)\n",
    "    ax_histx.plot(tmp_distribution_x,label=k,color=j)\n",
    "    ax_histy.plot(tmp_distribution_y,range(len(tmp_distribution_y)),label=k,color=j)\n",
    "ax.axis('off')\n",
    "#ax.set_title('Neocortex projection distribution',fontsize=fontsize_title*0.6)\n",
    "ax_histx.axis('off')\n",
    "ax_histy.axis('off')\n",
    "ax_histx.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#figure 5\n",
    "IC_a_m = -1000*1.21+5472\n",
    "IC_m_p = -1000*0.01+5472\n",
    "all_neuron_name = []\n",
    "for i in cortex_neuron:\n",
    "    all_neuron_name.append(i)\n",
    "xiao_region = ['ICa','ICm','ICp']\n",
    "xiao_region_color = {'ICa':'#DB4537','ICm':'#55AF7B','ICp':'#4286F3'}\n",
    "region_transform = {}\n",
    "for i in all_neuron_name:\n",
    "    if cortex_neuron[i]['soma_position'][0]<IC_a_m:\n",
    "        region = 'ICa'\n",
    "    else:\n",
    "        if cortex_neuron[i]['soma_position'][0]<IC_m_p:\n",
    "            region = 'ICm'\n",
    "        else:\n",
    "            region = 'ICp'\n",
    "    region_transform[i] = region\n",
    "def generate_statistics_information(neuron_list,region):\n",
    "    projection_information = {}\n",
    "    for i in tool_bar(neuron_list,desc='generating neuron information...',leave=False):\n",
    "        projection_information[i] = {}\n",
    "        projection_information[i]['axonal_projection_length_combined'] = get_axonal_statics_by_regions(i,region,'projectregion')\n",
    "        projection_information[i]['axonal_projection_length_left'] = get_axonal_statics_by_regions(i,region,'projectleftregion')\n",
    "        projection_information[i]['axonal_projection_length_right'] = get_axonal_statics_by_regions(i,region,'projectrightregion')\n",
    "        projection_information[i]['axonal_projection_terminal_combined'] = get_axonal_statics_by_regions(i,region,'terminalregion')\n",
    "        projection_information[i]['axonal_projection_terminal_left'] = get_axonal_statics_by_regions(i,region,'terminalleftregion')\n",
    "        projection_information[i]['axonal_projection_terminal_right'] = get_axonal_statics_by_regions(i,region,'terminalrightregion')\n",
    "    return projection_information"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#soma distribution\n",
    "flatmap = cv.imread('data/flatmapedge2.tif',1)\n",
    "fig = plt.figure(figsize=(9,4))\n",
    "ax = fig.add_subplot(1,2,1)\n",
    "group = {}\n",
    "for i in cortex_neuron:\n",
    "    region = region_transform[i]\n",
    "    position = cortex_neuron[i]['soma_flatten_position']\n",
    "    if region not in group:\n",
    "        group[region] = []\n",
    "    group[region].append(position)\n",
    "ax.matshow(flatmap)\n",
    "for i in group:\n",
    "    tmp = np.array(group[i])*2\n",
    "    ax.scatter(tmp[:,0],tmp[:,1],marker='o',color=xiao_region_color[i],s=0.5,label=i)\n",
    "ax.set_xticks([],[])\n",
    "ax.set_yticks([],[])\n",
    "ax.set_xlabel('L-M',fontsize=fontsize_label*0.6)\n",
    "ax.set_ylabel('P-A',fontsize=fontsize_label*0.6,rotation=90)\n",
    "ax.legend()\n",
    "ax = fig.add_subplot(1,2,2)\n",
    "group = {}\n",
    "for i in cortex_neuron:\n",
    "    region = region_transform[i]\n",
    "    position = cortex_neuron[i]['soma_position']\n",
    "    if region not in group:\n",
    "        group[region] = []\n",
    "    group[region].append(position)\n",
    "mu = 100  # mean of distribution\n",
    "sigma = 15  # standard deviation of distribution\n",
    "for i in group:\n",
    "    tmp = np.array(group[i])\n",
    "    ax.hist(tmp[:,0],density=True,color=xiao_region_color[i],label=i)\n",
    "ax.set_xticks([],[])\n",
    "ax.set_yticks([],[])\n",
    "ax.set_xlabel('A-P',fontsize=fontsize_label*0.6,rotation=0)\n",
    "ax.set_ylabel('Relative distribution',fontsize=fontsize_label*0.6,rotation=90)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Isocortex\n",
    "target_region_name = 'Isocortex'\n",
    "target_region_mask = get_region_mask(target_region_name)\n",
    "unit = 1\n",
    "sampling = 20\n",
    "kernel_size = 0\n",
    "projection_matrix_size = np.array(target_region_mask.shape).astype(int)//unit\n",
    "projection_matrix = {}\n",
    "for i in xiao_region:\n",
    "    projection_matrix[i] = {}\n",
    "    projection_matrix[i]['matrix'] = np.zeros([projection_matrix_size[0],projection_matrix_size[1]])\n",
    "    projection_matrix[i]['num'] = 0\n",
    "for i in tool_bar(all_neuron_name):\n",
    "    if cortex_neuron[i]['axonal_projection_length_combined'][:len(target_region_hierarchy['Cortex'])]==0:\n",
    "        continue\n",
    "    region = region_transform[i]\n",
    "    control_console_output(0)\n",
    "    neuron_tree = iondata.getNeuronTreeByID(i[:-7],i[-7:])\n",
    "    for point in neuron_tree.points[::sampling]:\n",
    "        if point.type==3:\n",
    "            continue\n",
    "        tmp = get_mirrored_point(point.xyz)\n",
    "        if if_inside_mask(target_region_mask,tmp/10):\n",
    "            a = np.array([np.round(tmp[0]/10/unit),np.round(tmp[1]/10/unit)]).astype(int)\n",
    "            b = np.where(a-kernel_size<0,0,a-kernel_size)\n",
    "            c = a+kernel_size+1\n",
    "            projection_matrix[region]['matrix'][b[0]:c[0],b[1]:c[1]] += 1\n",
    "    projection_matrix[region]['num'] += 1\n",
    "    control_console_output(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in xiao_region:\n",
    "    fig = plt.figure(figsize=(8, 5))\n",
    "    gs = fig.add_gridspec(2, 2,  width_ratios=(4.25, 0.75), height_ratios=(0.75, 3.25),left=0.1, right=0.9, bottom=0.1, top=0.9,wspace=0.05, hspace=0.05)\n",
    "    ax = fig.add_subplot(gs[1, 0])\n",
    "    ax_histx = fig.add_subplot(gs[0, 0], sharex=ax)\n",
    "    ax_histy = fig.add_subplot(gs[1, 1], sharey=ax)\n",
    "    tmp = projection_matrix[i]['matrix']/projection_matrix[i]['num']\n",
    "    tmp = np.log(tmp+1)\n",
    "    tmp /= np.max(tmp)\n",
    "    tmp += np.where(np.sum(target_region_mask,axis=-1)>0,0.02,0)\n",
    "    ax.matshow(tmp.T[0:570,150:1100],cmap='hot')\n",
    "    tmp_distribution = projection_matrix[i]['matrix'][150:1100,0:570]\n",
    "    tmp_distribution_x = np.sum(tmp_distribution,axis=1)\n",
    "    try:\n",
    "        tmp_distribution_x /= np.linalg.norm(tmp_distribution_x)\n",
    "    except:\n",
    "        pass\n",
    "    tmp_distribution_y = np.sum(tmp_distribution,axis=0)\n",
    "    try:\n",
    "        tmp_distribution_y /= np.linalg.norm(tmp_distribution_y)\n",
    "    except:\n",
    "        pass\n",
    "    ax_histx.plot(tmp_distribution_x,label=i,color=xiao_region_color[i])\n",
    "    ax_histy.plot(tmp_distribution_y,range(len(tmp_distribution_y)),label=i,color=xiao_region_color[i])\n",
    "    ax.axis('off')\n",
    "    #ax.set_title('Neocortex projection distribution',fontsize=fontsize_title*0.6)\n",
    "    ax_histx.axis('off')\n",
    "    ax_histy.axis('off')\n",
    "    ax_histx.legend()\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#amygdala\n",
    "_region_hierarchy = {\n",
    "        'OLF':['COA'],\n",
    "        'Cortex-sp':['LA','BLA','BMA','PA'],\n",
    "        'CNU':['AAA','BA','CEA','IA','MEA'],\n",
    "    }\n",
    "_region = []\n",
    "for i in _region_hierarchy:\n",
    "    for j in _region_hierarchy[i]:\n",
    "        _region.append(j)\n",
    "_projection_information = generate_statistics_information(all_neuron_name,_region)\n",
    "_projection_matrix = {}\n",
    "for i in _projection_information:\n",
    "    region = region_transform[i]\n",
    "    if region not in _projection_matrix:\n",
    "        _projection_matrix[region] = []\n",
    "    _projection_matrix[region].append(_projection_information[i]['axonal_projection_length_combined'])\n",
    "fig = plt.figure(figsize=(len(_region),4))\n",
    "ax = fig.add_subplot(1,1,1)\n",
    "for i,j in zip(range(3),['AIa','AIp','GI']):\n",
    "    #ax.bar([k*2+i/2 for k in range(len(_region))],np.log(np.mean(_projection_matrix[j],axis=0)+1),width=0.45,label=j,color=tmp_color[j],capstyle='round')\n",
    "    ax.bar([k*2+i/2 for k in range(len(_region))],np.power(np.mean(_projection_matrix[j],axis=0),0.5),width=0.45,label=j,color=tmp_color[j],capstyle='round')\n",
    "ax.set_xticks([i*2+0.5 for i in range(len(_region))],_region,fontsize=fontsize_ticks*0.6)\n",
    "#ax.set_yticks([0,50,100],[0,50,100],fontsize=fontsize_ticks*0.6)\n",
    "ax.set_xlabel('Target region',fontsize=fontsize_label*0.6)\n",
    "ax.set_ylabel('Axonal length ($\\sqrt{\\u00B5m}$)',fontsize=fontsize_label*0.6,rotation=90)\n",
    "ax.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#PAL\n",
    "_region_hierarchy = {\n",
    "        'PAL':['PALd','PALv','PALm','PALc'],\n",
    "    }\n",
    "_region = []\n",
    "for i in _region_hierarchy:\n",
    "    for j in _region_hierarchy[i]:\n",
    "        _region.append(j)\n",
    "_projection_information = generate_statistics_information(all_neuron_name,_region)\n",
    "_projection_matrix = {}\n",
    "for i in _projection_information:\n",
    "    region = region_transform[i]\n",
    "    if region not in _projection_matrix:\n",
    "        _projection_matrix[region] = []\n",
    "    _projection_matrix[region].append(_projection_information[i]['axonal_projection_length_combined'])\n",
    "fig = plt.figure(figsize=(len(_region),4))\n",
    "ax = fig.add_subplot(1,1,1)\n",
    "for i,j in zip(range(3),['AIa','AIp','GI']):\n",
    "    #ax.bar([k*2+i/2 for k in range(len(_region))],np.log(np.mean(_projection_matrix[j],axis=0)+1),width=0.45,label=j,color=tmp_color[j],capstyle='round')\n",
    "    ax.bar([k*2+i/2 for k in range(len(_region))],np.power(np.mean(_projection_matrix[j],axis=0),0.5),width=0.45,label=j,color=tmp_color[j],capstyle='round')\n",
    "ax.set_xticks([i*2+0.5 for i in range(len(_region))],_region,fontsize=fontsize_ticks*0.6)\n",
    "#ax.set_yticks([0,50,100],[0,50,100],fontsize=fontsize_ticks*0.6)\n",
    "ax.set_xlabel('Target region',fontsize=fontsize_label*0.6)\n",
    "ax.set_ylabel('Axonal length ($\\sqrt{\\u00B5m}$)',fontsize=fontsize_label*0.6,rotation=90)\n",
    "ax.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Striatum\n",
    "_region_hierarchy = {\n",
    "        'STRv':['ACB','FS','OT','LSS'],\n",
    "    }\n",
    "_region = []\n",
    "for i in _region_hierarchy:\n",
    "    for j in _region_hierarchy[i]:\n",
    "        _region.append(j)\n",
    "_projection_information = generate_statistics_information(all_neuron_name,_region)\n",
    "_projection_matrix = {}\n",
    "for i in _projection_information:\n",
    "    region = region_transform[i]\n",
    "    if region not in _projection_matrix:\n",
    "        _projection_matrix[region] = []\n",
    "    _projection_matrix[region].append(_projection_information[i]['axonal_projection_length_combined'])\n",
    "fig = plt.figure(figsize=(len(_region),4))\n",
    "ax = fig.add_subplot(1,1,1)\n",
    "for i,j in zip(range(3),['AIa','AIp','GI']):\n",
    "    #ax.bar([k*2+i/2 for k in range(len(_region))],np.log(np.mean(_projection_matrix[j],axis=0)+1),width=0.45,label=j,color=tmp_color[j],capstyle='round')\n",
    "    ax.bar([k*2+i/2 for k in range(len(_region))],np.power(np.mean(_projection_matrix[j],axis=0),0.5),width=0.45,label=j,color=tmp_color[j],capstyle='round')\n",
    "ax.set_xticks([i*2+0.5 for i in range(len(_region))],_region,fontsize=fontsize_ticks*0.6)\n",
    "#ax.set_yticks([0,50,100],[0,50,100],fontsize=fontsize_ticks*0.6)\n",
    "ax.set_xlabel('Target region',fontsize=fontsize_label*0.6)\n",
    "ax.set_ylabel('Axonal length ($\\sqrt{\\u00B5m}$)',fontsize=fontsize_label*0.6,rotation=90)\n",
    "ax.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#NTS\n",
    "nts_mask = get_region_mask('NTS')\n",
    "unit = 1\n",
    "sampling = 20\n",
    "kernel_size = 3\n",
    "projection_matrix_size = np.array(nts_mask.shape).astype(int)//unit\n",
    "projection_matrix = {}\n",
    "for i in ['AIa','AIp','GI']:\n",
    "    slected_neuron_region_index.append(i)\n",
    "    projection_matrix[i] = {}\n",
    "    projection_matrix[i]['matrix'] = np.zeros([projection_matrix_size[0],projection_matrix_size[2]])\n",
    "    projection_matrix[i]['num'] = 0\n",
    "for i in tool_bar(all_neuron_name):\n",
    "    if cortex_neuron[i]['axonal_projection_length_combined'][target_region.index('NTS')]==0:\n",
    "        continue\n",
    "    region = region_transform[i]\n",
    "    control_console_output(0)\n",
    "    neuron_tree = iondata.getNeuronTreeByID(i[:-7],i[-7:])\n",
    "    for point in neuron_tree.points[::sampling]:\n",
    "        if point.type==3:\n",
    "            continue\n",
    "        tmp = get_mirrored_point(point.xyz)\n",
    "        if if_inside_mask(nts_mask,tmp/10):\n",
    "            a = np.array([np.round(tmp[0]/10/unit),np.round(tmp[2]/10/unit)]).astype(int)\n",
    "            b = np.where(a-kernel_size<0,0,a-kernel_size)\n",
    "            c = a+kernel_size+1\n",
    "            projection_matrix[region]['matrix'][b[0]:c[0],b[1]:c[1]] += 1\n",
    "    projection_matrix[region]['num'] += 1\n",
    "    control_console_output(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(10, 8))\n",
    "gs = fig.add_gridspec(2, 2,  width_ratios=(4.25, 0.75), height_ratios=(0.75, 3.25),left=0.1, right=0.9, bottom=0.1, top=0.9,wspace=0.05, hspace=0.05)\n",
    "ax = fig.add_subplot(gs[1, 0])\n",
    "ax_histx = fig.add_subplot(gs[0, 0], sharex=ax)\n",
    "ax_histy = fig.add_subplot(gs[1, 1], sharey=ax)\n",
    "tmp = []\n",
    "for i in ['AIa','GI','AIp']:\n",
    "    tmp.append(projection_matrix[i]['matrix'])\n",
    "    tmp[-1] = np.log(tmp[-1]+1)\n",
    "    #tmp[-1] = np.power(tmp[-1],0.5)\n",
    "    tmp[-1] /= np.max(tmp[-1])\n",
    "    #tmp[-1] = np.where(tmp[-1]<0.9,0,tmp[-1])\n",
    "tmp = np.array(tmp).transpose(1,2,0)\n",
    "ax.matshow(tmp[1000:,360:760])\n",
    "for i,j,k in zip(projection_matrix,['r','b','g'],['AIa','AIp','GI']):\n",
    "    tmp_distribution = projection_matrix[i]['matrix'][1000:,360:760]\n",
    "    tmp_distribution_x = np.sum(tmp_distribution,axis=0)\n",
    "    try:\n",
    "        tmp_distribution_x /= np.linalg.norm(tmp_distribution_x)\n",
    "    except:\n",
    "        pass\n",
    "    tmp_distribution_y = np.sum(tmp_distribution,axis=1)\n",
    "    try:\n",
    "        tmp_distribution_y /= np.linalg.norm(tmp_distribution_y)\n",
    "    except:\n",
    "        pass\n",
    "    ax_histx.plot(tmp_distribution_x,label=k,color=j)\n",
    "    ax_histy.plot(tmp_distribution_y,range(len(tmp_distribution_y)),label=k,color=j)\n",
    "ax.axis('off')\n",
    "#ax.set_title('Neocortex projection distribution',fontsize=fontsize_title*0.6)\n",
    "ax_histx.axis('off')\n",
    "ax_histy.axis('off')\n",
    "ax_histx.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in ['AIa','AIp','GI']:\n",
    "    plt.imshow(projection_matrix[i]['matrix'][1000:,350:750])\n",
    "    plt.colorbar()\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# brain netwworks (allen)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_cortex_neuron_list = []\n",
    "samples = iondata.getSampleInfo(projectID='ISO')\n",
    "for sample in samples:\n",
    "    for neuron in iondata.getNeuronListBySampleID(str(sample['fMOST_id'])):\n",
    "        all_cortex_neuron_list.append(neuron['sampleid']+neuron['name'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tmp_transform_1,tmp_transform_2 = get_subregions(target_region_hierarchy['Cortex'])\n",
    "core_region = ['AId','AIv','AIp','GU','VISC']\n",
    "tmplate = []\n",
    "for i in core_region:\n",
    "    tmplate.append(0.0)\n",
    "tmplate = np.array(tmplate)\n",
    "IC_mask = get_region_mask(core_region[0])\n",
    "for i,j in enumerate(core_region[1:]):\n",
    "    IC_mask = np.where(get_region_mask(j),i+2,IC_mask)\n",
    "insula_input = {}\n",
    "insula_input_count = {}\n",
    "insula_input_terminal = {}\n",
    "for neuron in tool_bar(all_cortex_neuron_list):\n",
    "    control_console_output(0)\n",
    "    neuron_tree = iondata.getNeuronTreeByID(neuron[:-7],neuron[-7:])\n",
    "    neuron_property = iondata.getNeuronPropertyByID(neuron[:-7],neuron[-7:])\n",
    "    control_console_output(1)\n",
    "    try:\n",
    "        region = tmp_transform_2[neuron_property['somaregion']]\n",
    "    except:\n",
    "        continue\n",
    "    if region in ['AI']:\n",
    "        try:\n",
    "            region += cortex_neuron[neuron]['cortex_subregion']\n",
    "        except:\n",
    "            continue\n",
    "    if region not in insula_input:\n",
    "        insula_input[region] = tmplate*0\n",
    "        insula_input_count[region] = 0\n",
    "        insula_input_terminal[region] = tmplate*0\n",
    "    for edge in neuron_tree.edges:\n",
    "        if edge.data[0].type==3:\n",
    "            continue\n",
    "        start = np.array(edge.data[0].xyz)\n",
    "        for point in edge.data[1:]:\n",
    "            end = np.array(point.xyz)\n",
    "            value = int(get_mask_value(IC_mask,end/10))\n",
    "            if value>0:\n",
    "                length = np.linalg.norm(end-start)\n",
    "                insula_input[region][value-1] += length\n",
    "            start = end\n",
    "    for terminal in neuron_tree.terminals:\n",
    "        if terminal.type==3:\n",
    "            continue\n",
    "        point = np.array(terminal.xyz)\n",
    "        value = int(get_mask_value(IC_mask,point/10))\n",
    "        if value>0:\n",
    "            insula_input_terminal[region][value-1] += 1\n",
    "    insula_input_count[region] += 1\n",
    "for i in insula_input:\n",
    "    insula_input[i] /= insula_input_count[i]\n",
    "    insula_input_terminal[i] /= insula_input_count[i]\n",
    "insula_output = {}\n",
    "insula_output_count = tmplate*0\n",
    "insula_output_terminal = {}\n",
    "for neuron in cortex_neuron:\n",
    "    region = cortex_neuron[neuron]['cortex_region']\n",
    "    if region in ['AI']:\n",
    "        region += cortex_neuron[neuron]['cortex_subregion']\n",
    "    for i,j in zip(target_region,cortex_neuron[neuron]['axonal_projection_length_combined']):\n",
    "        if i not in insula_output:\n",
    "            insula_output[i] = tmplate*0\n",
    "            insula_output_terminal[i] = tmplate*0\n",
    "        insula_output[i][core_region.index(region)] += j\n",
    "    insula_output_count[core_region.index(region)] += 1\n",
    "for neuron in cortex_neuron:\n",
    "    region = cortex_neuron[neuron]['cortex_region']\n",
    "    if region in ['AI']:\n",
    "        region += cortex_neuron[neuron]['cortex_subregion']\n",
    "    for i,j in zip(target_region,cortex_neuron[neuron]['axonal_projection_terminal_combined']):\n",
    "        insula_output_terminal[i][core_region.index(region)] += j\n",
    "for i in insula_output:\n",
    "    insula_output[i] /= insula_output_count\n",
    "    insula_output_terminal[i] /= insula_output_count"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.save(data_path+'insula_input_allen.npy',insula_input)\n",
    "np.save(data_path+'insula_input_terminal_allen.npy',insula_input_terminal)\n",
    "np.save(data_path+'insula_output_allen.npy',insula_output)\n",
    "np.save(data_path+'insula_output_terminal_allen.npy',insula_output_terminal)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "insula_input = np.load(data_path+'insula_input_allen.npy',allow_pickle=True).item()\n",
    "insula_input_terminal = np.load(data_path+'insula_input_terminal_allen.npy',allow_pickle=True).item()\n",
    "insula_output = np.load(data_path+'insula_output_allen.npy',allow_pickle=True).item()\n",
    "insula_output_terminal = np.load(data_path+'insula_output_terminal_allen.npy',allow_pickle=True).item()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# unnormalized version"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tmp_transform_1,tmp_transform_2 = get_subregions(target_region_hierarchy['Cortex'])\n",
    "core_region = ['AId','AIv','AIp','GU','VISC']\n",
    "tmplate = []\n",
    "for i in core_region:\n",
    "    tmplate.append(0.0)\n",
    "tmplate = np.array(tmplate)\n",
    "insula_input_count = {}\n",
    "for neuron in tool_bar(all_cortex_neuron_list):\n",
    "    control_console_output(0)\n",
    "    neuron_property = iondata.getNeuronPropertyByID(neuron[:-7],neuron[-7:])\n",
    "    control_console_output(1)\n",
    "    try:\n",
    "        region = tmp_transform_2[neuron_property['somaregion']]\n",
    "    except:\n",
    "        continue\n",
    "    if region in ['AI']:\n",
    "        try:\n",
    "            region += cortex_neuron[neuron]['cortex_subregion']\n",
    "        except:\n",
    "            continue\n",
    "    if region not in insula_input_count:\n",
    "        insula_input_count[region] = 0\n",
    "    insula_input_count[region] += 1\n",
    "insula_output_count = tmplate*0\n",
    "for neuron in cortex_neuron:\n",
    "    region = cortex_neuron[neuron]['cortex_region']\n",
    "    if region in ['AI']:\n",
    "        region += cortex_neuron[neuron]['cortex_subregion']\n",
    "    insula_output_count[core_region.index(region)] += 1\n",
    "for i in insula_input:\n",
    "    insula_input[i] *= insula_input_count[i]\n",
    "    insula_input_terminal[i] *= insula_input_count[i]\n",
    "for i in insula_output:\n",
    "    insula_output[i] *= insula_output_count\n",
    "    insula_output_terminal[i] *= insula_output_count"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "_insula_input = {}\n",
    "for i in insula_input:\n",
    "    _insula_input[i] = insula_input[i][:4]\n",
    "    _insula_input[i][3] += insula_input[i][4]\n",
    "insula_input = _insula_input\n",
    "insula_input['GI'] = (insula_input['GU']+insula_input['VISC'])/2\n",
    "del insula_input['GU']\n",
    "del insula_input['VISC']\n",
    "_insula_input_terminal = {}\n",
    "for i in insula_input_terminal:\n",
    "    _insula_input_terminal[i] = insula_input_terminal[i][:4]\n",
    "    _insula_input_terminal[i][3] += insula_input_terminal[i][4]\n",
    "insula_input_terminal = _insula_input_terminal\n",
    "insula_input_terminal['GI'] = (insula_input_terminal['GU']+insula_input_terminal['VISC'])/2\n",
    "del insula_input_terminal['GU']\n",
    "del insula_input_terminal['VISC']\n",
    "_insula_output = {}\n",
    "for i in insula_output:\n",
    "    _insula_output[i] = insula_output[i][:4]\n",
    "    _insula_output[i][3] += insula_output[i][4]\n",
    "insula_output = _insula_output\n",
    "insula_output['GI'] = (insula_output['GU']+insula_output['VISC'])/2\n",
    "del insula_output['GU']\n",
    "del insula_output['VISC']\n",
    "del insula_output['AI']\n",
    "_insula_output_terminal = {}\n",
    "for i in insula_output_terminal:\n",
    "    _insula_output_terminal[i] = insula_output_terminal[i][:4]\n",
    "    _insula_output_terminal[i][3] += insula_output_terminal[i][4]\n",
    "insula_output_terminal = _insula_output_terminal\n",
    "insula_output_terminal['GI'] = (insula_output_terminal['GU']+insula_output_terminal['VISC'])/2\n",
    "del insula_output_terminal['GU']\n",
    "del insula_output_terminal['VISC']\n",
    "del insula_output_terminal['AI']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.save(data_path+'insula_input.npy',insula_input)\n",
    "np.save(data_path+'insula_output.npy',insula_output)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import networkx as nx\n",
    "from netgraph import Graph\n",
    "threshold = 100000*3.5\n",
    "node_layout = {\n",
    "    'ORB':np.array([-4.25-1,0.25]),'ILA':np.array([-3.5-1,1.75]),'PL':np.array([-3.5-1,1.25]),'FRP':np.array([-3.75-1,1]),'ACA':np.array([-3.0-1,1.75]),'RSP':np.array([-2.5-1,1.75]),'PTLp':np.array([-2.0-1,1.75]),'VIS':np.array([-1.5-1,1.75]),\n",
    "    'MOs':np.array([-2-1,2]),'MOp':np.array([-1.5-1,2]),'SSp':np.array([-1.0-1,2]),'SSs':np.array([-0.5-1,2]),\n",
    "    'AUD':np.array([-0.0-1,2]),'TEa':np.array([-0.5-1,2]),'ECT':np.array([1-1,2]),'PERI':np.array([1.5-1,2]),\n",
    "    'AId':np.array([-0.4,-0.1]),'AIv':np.array([-0.4,-0.4]),'AIp':np.array([0.4,-0.3]),'GI':np.array([0.0,0.2]),\n",
    "    'AON':np.array([-3.25,-2]),'PIR':np.array([-3,-2]),'TR':np.array([-2.5,-2.0]),'NLOT':np.array([-2.25,-1.5]),'PAA':np.array([-2,-1.5]),'CTXsp':np.array([-2,-2]),'RHP':np.array([-1,-1.5]),\n",
    "    'PAL':np.array([0.25,-2]),'STRd':np.array([1.5,-1.5]),'STRv':np.array([1,-2]),'sAMY':np.array([-0.5,-2]),'COA':np.array([-1,-2]),\n",
    "    'MED':np.array([2,2]),'MTN':np.array([2,1.75]),'ILM':np.array([2,1.5]),'VENT':np.array([2,1]),'LAT':np.array([2,0.5]),\n",
    "    'LZ':np.array([3,-0.2]),\n",
    "    'SNr':np.array([3.75,-0.75]),'PAG':np.array([4,-0.5]),'MRN':np.array([4,-1]),\n",
    "    'PG':np.array([6,-0.5]),'PB':np.array([6,0.5]),\n",
    "    'PARN':np.array([6,0]),'NTS':np.array([6.5,0.25]),\n",
    "}\n",
    "brainstem_region_volume_center_range = np.load(data_path+'brainstem_region_volume_center_range.npy',allow_pickle=True).item()\n",
    "node_layout = {}\n",
    "for i in brainstem_region_volume_center_range:\n",
    "    node_layout[i] = np.array(brainstem_region_volume_center_range[i][1][:2])/100\n",
    "    node_layout[i][0] *= 2\n",
    "    node_layout[i][1] = -node_layout[i][1]\n",
    "node_layout['AId'] = np.array([8.58724983-0.3, -4.66861863+0.3])\n",
    "node_layout['AIv'] = np.array([8.58724983-0.3, -4.66861863-0.3])\n",
    "node_layout['AIp'] = np.array([9.58724983, -4.66861863])\n",
    "node_layout['GI'] = np.array([8.9612381 , -4.45964817])\n",
    "node_layout['LZ'] = node_layout['ZI']\n",
    "node_layout['CTXsp'] = node_layout['EP']\n",
    "node_layout['PVR'] = node_layout['SBPV']\n",
    "node_layout['PVZ'] = node_layout['PVi']\n",
    "node_layout['MEZ'] = node_layout['AHN']\n",
    "node_layout['STRv'] = node_layout['OT']\n",
    "node_layout['PAL'] = node_layout['PALd']\n",
    "scale = threshold*25/3.5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import cv2\n",
    "root_mask = get_region_mask('root')\n",
    "root_mask = np.sum(root_mask,axis=-1)\n",
    "root_mask = np.where(root_mask>0,1,0)\n",
    "new_width = root_mask.shape[1]\n",
    "new_height = root_mask.shape[0]*2\n",
    "root_mask = cv2.resize(root_mask.astype('uint8'),(new_width, new_height))\n",
    "contours, hierarchy = cv2.findContours(root_mask, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE)\n",
    "root_mask = cv2.drawContours(root_mask*0, contours, -1,1)\n",
    "plt.imshow(1-root_mask.T,cmap='gray')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## input"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mode = 0\n",
    "color_mode = 1\n",
    "insula_network = nx.DiGraph()\n",
    "core_region = []\n",
    "node_color = {}\n",
    "node_color2 = {'AId':np.array([0.859,0.267,0.216]),'AIv':np.array([0.059,0.616,0.345]),'AIp':np.array([0.258,0.52,0.957]),'GI':np.array([0.957,0.627,0.0])}\n",
    "node_labels = {}\n",
    "node_size = {}\n",
    "node_name = []\n",
    "for i in node_color2:\n",
    "    node_name.append(i)\n",
    "    core_region.append(i)\n",
    "    node_labels[i] = i\n",
    "    node_size[i] = 5\n",
    "    node_color[i] = rgb_to_hex(np.round(node_color2[i]*255).astype(int))\n",
    "node_color_count = {}\n",
    "edge_color = {}\n",
    "edge_width = {}\n",
    "for i in node_name:\n",
    "    insula_network.add_node(i)\n",
    "for i in insula_input:\n",
    "    if np.sum(insula_input[i])<threshold:\n",
    "        continue\n",
    "    if i in core_region:\n",
    "        continue\n",
    "    tmp = 0\n",
    "    for j,k in zip(core_region,insula_input[i]):\n",
    "        if i not in node_name:\n",
    "            node_name.append(i)\n",
    "            insula_network.add_node(i)\n",
    "            if color_mode==0:\n",
    "                node_color[i] = '#fac230'\n",
    "            else:\n",
    "                node_color[i] = np.array([0.0 for _i in core_region])\n",
    "        insula_network.add_edge(i,j)\n",
    "        if mode==0:\n",
    "            edge_width[i,j] = k/scale\n",
    "        else:\n",
    "            edge_width[i,j] = np.log(k)/3\n",
    "        tmp += edge_width[i,j]\n",
    "        edge_color[i,j] = node_color[j]\n",
    "        if color_mode==1:\n",
    "            node_color[i][core_region.index(j)] += edge_width[i,j]\n",
    "    node_size[i] = tmp**0.5*5\n",
    "if color_mode==1:\n",
    "    for i in node_color:\n",
    "        if i in core_region:\n",
    "            continue\n",
    "        factors = node_color[i]/np.sum(node_color[i])\n",
    "        node_color[i] =np.sum([_i*node_color2[_j] for _i,_j in zip(factors,core_region)],axis=0)\n",
    "        node_color[i] = rgb_to_hex(np.round(node_color[i]*255).astype(int))\n",
    "for i in insula_input:\n",
    "    if i not in ['ICa','ICm','ICp']:\n",
    "        continue\n",
    "    for j,k in zip(['ICa','ICm','ICp'],insula_input[i]):\n",
    "        insula_network.add_edge(i,j)\n",
    "        if mode==0:\n",
    "            edge_width[i,j] = k/1000\n",
    "        else:\n",
    "            edge_width[i,j] = np.log(k)/3\n",
    "        edge_color[i,j] = node_color[i]\n",
    "fig = plt.figure(figsize=(12,12))\n",
    "ax = plt.subplot(1,1,1)\n",
    "ax.imshow(1-root_mask.T,extent=(0,root_mask.shape[0]/100,-root_mask.shape[1]/100,0),cmap='gray')\n",
    "tmp = Graph(insula_network,\n",
    "      node_size=node_size,\n",
    "      node_color=node_color, \n",
    "      node_edge_width=1,\n",
    "      node_labels=True,\n",
    "      node_label_fontdict=dict(size=7),\n",
    "      node_label_offset=0.0,\n",
    "      edge_color=edge_color,  \n",
    "      edge_width=edge_width,        \n",
    "      edge_alpha=0.75,\n",
    "      node_layout = node_layout,\n",
    "      edge_layout='bundled',\n",
    "      edge_layout_kwargs={'k':250,'total_cycles':8,'total_iterations':50,'compatibility_threshold':0.0025},\n",
    "      with_labels=True,\n",
    "      ax=ax,\n",
    ")\n",
    "#ax.set_xlim(-6,7)\n",
    "#ax.set_ylim(-3,3)\n",
    "plt.savefig(data_path+'IC_input_allen_u.svg',format='svg')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mode = 0\n",
    "color_mode = 1\n",
    "insula_network = nx.DiGraph()\n",
    "core_region = []\n",
    "node_color = {}\n",
    "node_color2 = {'AId':np.array([0.859,0.267,0.216]),'AIv':np.array([0.059,0.616,0.345]),'AIp':np.array([0.258,0.52,0.957]),'GI':np.array([0.957,0.627,0.0])}\n",
    "node_labels = {}\n",
    "node_size = {}\n",
    "node_name = []\n",
    "for i in node_color2:\n",
    "    node_name.append(i)\n",
    "    core_region.append(i)\n",
    "    node_labels[i] = i\n",
    "    node_size[i] = 5\n",
    "    node_color[i] = rgb_to_hex(np.round(node_color2[i]*255).astype(int))\n",
    "node_color_count = {}\n",
    "edge_color = {}\n",
    "edge_width = {}\n",
    "for i in node_name:\n",
    "    insula_network.add_node(i)\n",
    "for i in insula_input:\n",
    "    if i not in core_region:\n",
    "        continue\n",
    "    for j,k in zip(core_region,insula_input[i]):\n",
    "        insula_network.add_edge(i,j)\n",
    "        if mode==0:\n",
    "            edge_width[i,j] = k/scale\n",
    "        else:\n",
    "            edge_width[i,j] = np.log(k)/3\n",
    "        edge_color[i,j] = node_color[i]\n",
    "for i in insula_output:\n",
    "    if np.sum(insula_output[i])<threshold:\n",
    "        continue\n",
    "    if i in core_region:\n",
    "        continue\n",
    "    tmp = 0\n",
    "    for j,k in zip(core_region,insula_output[i]):\n",
    "        if i not in node_name:\n",
    "            node_name.append(i)\n",
    "            if color_mode==0:\n",
    "                node_color[i] = '#f032e6'\n",
    "            else:\n",
    "                node_color[i] = np.array([0.0 for _i in core_region])\n",
    "            insula_network.add_node(i)\n",
    "        insula_network.add_edge(j,i)\n",
    "        if mode==0:\n",
    "            edge_width[j,i] = k/scale\n",
    "            if i=='STRd':\n",
    "                edge_width[j,i] /= 1\n",
    "        else:\n",
    "            edge_width[j,i] = np.log(k)/3\n",
    "        tmp += edge_width[j,i]\n",
    "        edge_color[j,i] = node_color[j]\n",
    "        if color_mode==1:\n",
    "            node_color[i][core_region.index(j)] += edge_width[j,i]\n",
    "    node_size[i] = tmp**0.5*5\n",
    "if color_mode==1:\n",
    "    for i in node_color:\n",
    "        if i in core_region:\n",
    "            continue\n",
    "        factors = node_color[i]/np.sum(node_color[i])\n",
    "        node_color[i] =np.sum([_i*node_color2[_j] for _i,_j in zip(factors,core_region)],axis=0)\n",
    "        node_color[i] = rgb_to_hex(np.round(node_color[i]*255).astype(int))\n",
    "fig = plt.figure(figsize=(12,12))\n",
    "ax = plt.subplot(1,1,1)\n",
    "ax.imshow(1-root_mask.T,extent=(0,root_mask.shape[0]/100,-root_mask.shape[1]/100,0),cmap='gray')\n",
    "tmp = Graph(insula_network,\n",
    "      node_size=node_size,\n",
    "      node_color=node_color, \n",
    "      node_edge_width=1,\n",
    "      node_labels=True,\n",
    "      node_label_fontdict=dict(size=7),\n",
    "      node_label_offset=0.0,\n",
    "      edge_color=edge_color,  \n",
    "      edge_width=edge_width,        \n",
    "      edge_alpha=0.75,\n",
    "      node_layout = node_layout,\n",
    "      edge_layout='bundled',\n",
    "      edge_layout_kwargs={'k':250,'total_cycles':8,'total_iterations':50,'compatibility_threshold':0.0025},\n",
    "      with_labels=True,\n",
    "      ax=ax,\n",
    ")\n",
    "#ax.set_xlim(-6,7)\n",
    "#ax.set_ylim(-3,3)\n",
    "plt.savefig(data_path+'IC_ioutput_allen_u.svg',format='svg')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mode = 0\n",
    "color_mode = 1\n",
    "insula_network = nx.DiGraph()\n",
    "core_region = []\n",
    "node_color = {}\n",
    "node_color2 = {'AId':np.array([0.859,0.267,0.216]),'AIv':np.array([0.059,0.616,0.345]),'AIp':np.array([0.258,0.52,0.957]),'GI':np.array([0.957,0.627,0.0])}\n",
    "node_labels = {}\n",
    "node_size = {}\n",
    "node_name = []\n",
    "for i in node_color2:\n",
    "    node_name.append(i)\n",
    "    core_region.append(i)\n",
    "    node_labels[i] = i\n",
    "    node_size[i] = 5\n",
    "    node_color[i] = rgb_to_hex(np.round(node_color2[i]*255).astype(int))\n",
    "node_color_count = {}\n",
    "edge_color = {}\n",
    "edge_width = {}\n",
    "for i in node_name:\n",
    "    insula_network.add_node(i)\n",
    "input_node_tmp = []\n",
    "for i in insula_input:\n",
    "    tmp_strength = np.sum(insula_input[i])\n",
    "    if i in insula_output:\n",
    "        tmp_strength += np.sum(insula_output[i])\n",
    "    if tmp_strength<threshold:\n",
    "        continue\n",
    "    if i in core_region:\n",
    "        continue\n",
    "    tmp = 0\n",
    "    mark = 0\n",
    "    for j,k in zip(core_region,insula_input[i]):\n",
    "        if i not in node_name:\n",
    "            node_name.append(i)\n",
    "            insula_network.add_node(i)\n",
    "            if color_mode==0:\n",
    "                node_color[i] = '#fac230'\n",
    "            else:\n",
    "                node_color[i] = np.array([0.0 for _i in core_region])\n",
    "        insula_network.add_edge(i,j)\n",
    "        if mode==0:\n",
    "            edge_width[i,j] = k/scale\n",
    "        else:\n",
    "            edge_width[i,j] = np.log(k)/3\n",
    "        tmp += edge_width[i,j]\n",
    "        edge_color[i,j] = '#fac230'\n",
    "        if color_mode==1:\n",
    "            node_color[i][core_region.index(j)] += edge_width[i,j]\n",
    "        mark = 1\n",
    "    if mark:\n",
    "        input_node_tmp.append(i)\n",
    "        node_size[i] = tmp\n",
    "for i in insula_input:\n",
    "    if i not in core_region:\n",
    "        continue\n",
    "    for j,k in zip(core_region,insula_input[i]):\n",
    "        insula_network.add_edge(i,j)\n",
    "        if mode==0:\n",
    "            edge_width[i,j] = k/scale\n",
    "        else:\n",
    "            edge_width[i,j] = np.log(k)/3\n",
    "        edge_color[i,j] = node_color[i]\n",
    "for i in insula_output:\n",
    "    if np.sum(insula_output[i])<threshold:\n",
    "        continue\n",
    "    if i in core_region:\n",
    "        continue\n",
    "    tmp = 0\n",
    "    for j,k in zip(core_region,insula_output[i]):\n",
    "        if i not in node_name:\n",
    "            node_name.append(i)\n",
    "            insula_network.add_node(i)\n",
    "            if color_mode==0:\n",
    "                node_color[i] = '#f032e6'\n",
    "            else:\n",
    "                if i in input_node_tmp:\n",
    "                    node_color[i] += insula_output[i]\n",
    "                else:\n",
    "                    node_color[i] = np.array([0.0 for _i in core_region])\n",
    "        insula_network.add_edge(j,i)\n",
    "        if mode==0:\n",
    "            edge_width[j,i] = k/scale\n",
    "            if i=='STRd':\n",
    "                edge_width[j,i] /= 1\n",
    "        else:\n",
    "            edge_width[j,i] = np.log(k)/3\n",
    "        tmp += edge_width[j,i]\n",
    "        edge_color[j,i] = node_color[j]\n",
    "        if color_mode==1:\n",
    "            node_color[i][core_region.index(j)] += edge_width[j,i]\n",
    "    if i in input_node_tmp:\n",
    "        node_size[i] += tmp\n",
    "    else:\n",
    "        node_size[i] = tmp**0.5*5\n",
    "for i in input_node_tmp:\n",
    "    node_size[i] = node_size[i]**0.5*5\n",
    "    factors = node_color[i]/np.sum(node_color[i])\n",
    "if color_mode==1:\n",
    "    for i in node_color:\n",
    "        if i in core_region:\n",
    "            continue\n",
    "        factors = node_color[i]/np.sum(node_color[i])\n",
    "        node_color[i] =np.sum([_i*node_color2[_j] for _i,_j in zip(factors,core_region)],axis=0)\n",
    "        node_color[i] = rgb_to_hex(np.round(node_color[i]*255).astype(int))\n",
    "        for j in edge_color:\n",
    "            if j[0]==i:\n",
    "                edge_color[j] = node_color[i]\n",
    "fig = plt.figure(figsize=(12,12))\n",
    "ax = plt.subplot(1,1,1)\n",
    "ax.imshow(1-root_mask.T,extent=(0,root_mask.shape[0]/100,-root_mask.shape[1]/100,0),cmap='gray')\n",
    "tmp = Graph(insula_network,\n",
    "      node_size=node_size,\n",
    "      node_color=node_color, \n",
    "      node_edge_width=0.2,\n",
    "      node_labels=True,\n",
    "      node_label_fontdict=dict(size=7),\n",
    "      node_label_offset=0.0,\n",
    "      edge_color=edge_color,  \n",
    "      edge_width=edge_width,        \n",
    "      edge_alpha=0.75,        \n",
    "      #node_layout='spring',\n",
    "      #node_layout='circular',\n",
    "      #node_layout_kwargs=dict(subsets=node_sets),\n",
    "      node_layout = node_layout,\n",
    "      edge_layout='bundled',\n",
    "      edge_layout_kwargs={'k':250,'total_cycles':8,'total_iterations':50,'compatibility_threshold':0.0025},\n",
    "      with_labels=True,\n",
    "      ax=ax,\n",
    ")\n",
    "#ax.set_xlim(-6,7)\n",
    "#ax.set_ylim(-3,3)\n",
    "plt.savefig(data_path+'IC_inouput_allen_u.svg',format='svg')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "index = target_region.index('NTS')\n",
    "info = [0,0]\n",
    "for i in insular_neuron_name:\n",
    "    info[1] += cortex_neuron[i]['axonal_projection_length_combined'][index]\n",
    "    if cortex_neuron[i]['cortex_region']=='AI' and cortex_neuron[i]['cortex_subregion']=='d':\n",
    "        info[0] += cortex_neuron[i]['axonal_projection_length_combined'][index]\n",
    "print(info[0]/np.sum(info[1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# brain networks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_cortex_neuron_list = []\n",
    "samples = iondata.getSampleInfo(projectID='ISO')\n",
    "for sample in samples:\n",
    "    for neuron in iondata.getNeuronListBySampleID(str(sample['fMOST_id'])):\n",
    "        all_cortex_neuron_list.append(neuron['sampleid']+neuron['name'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tmp_transform_1,tmp_transform_2 = get_subregions(target_region_hierarchy['Cortex'])\n",
    "IC_mask = get_region_mask('AI')+get_region_mask('VISC')+get_region_mask('GU')\n",
    "IC_mask = np.where(IC_mask>0,1.0,0.0)\n",
    "insula_input = {}\n",
    "insula_input_count = {}\n",
    "insula_input_terminal = {}\n",
    "for neuron in tool_bar(all_cortex_neuron_list):\n",
    "    control_console_output(0)\n",
    "    neuron_tree = iondata.getNeuronTreeByID(neuron[:-7],neuron[-7:])\n",
    "    neuron_property = iondata.getNeuronPropertyByID(neuron[:-7],neuron[-7:])\n",
    "    control_console_output(1)\n",
    "    try:\n",
    "        region = tmp_transform_2[neuron_property['somaregion']]\n",
    "    except:\n",
    "        continue\n",
    "    if region in ['AI','GU','VISC']:\n",
    "        if neuron_property['somapoint'][0]<IC_a_m:\n",
    "            region = 'ICa'\n",
    "        else:\n",
    "            if neuron_property['somapoint'][0]<IC_m_p:\n",
    "                region = 'ICm'\n",
    "            else:\n",
    "                region = 'ICp'\n",
    "    if region not in insula_input:\n",
    "        insula_input[region] = np.array([0.0,0.0,0.0])\n",
    "        insula_input_count[region] = 0\n",
    "        insula_input_terminal[region] = np.array([0.0,0.0,0.0])\n",
    "    for edge in neuron_tree.edges:\n",
    "        if edge.data[0].type==3:\n",
    "            continue\n",
    "        start = np.array(edge.data[0].xyz)\n",
    "        for point in edge.data[1:]:\n",
    "            end = np.array(point.xyz)\n",
    "            if if_inside_mask(IC_mask,end/10):\n",
    "                length = np.linalg.norm(end-start)\n",
    "                if end[0]<IC_a_m:\n",
    "                    insula_input[region][0] += length\n",
    "                else:\n",
    "                    if end[0]<IC_m_p:\n",
    "                        insula_input[region][1] += length\n",
    "                    else:\n",
    "                        insula_input[region][2] += length\n",
    "            start = end\n",
    "    for terminal in neuron_tree.terminals:\n",
    "        if terminal.type==3:\n",
    "            continue\n",
    "        point = np.array(terminal.xyz)\n",
    "        if if_inside_mask(IC_mask,point/10):\n",
    "            if point[0]<IC_a_m:\n",
    "                insula_input_terminal[region][0] += 1\n",
    "            else:\n",
    "                if point[0]<IC_m_p:\n",
    "                    insula_input_terminal[region][1] += 1\n",
    "                else:\n",
    "                    insula_input_terminal[region][2] += 1\n",
    "    insula_input_count[region] += 1\n",
    "for i in insula_input:\n",
    "    insula_input[i] /= insula_input_count[i]\n",
    "    insula_input_terminal[i] /= insula_input_count[i]\n",
    "insula_output = {}\n",
    "insula_output_count = np.array([0,0,0])\n",
    "insula_output_terminal = {}\n",
    "for neuron in cortex_neuron:\n",
    "    for i,j in zip(target_region,cortex_neuron[neuron]['axonal_projection_length_combined']):\n",
    "        if i not in insula_output:\n",
    "            insula_output[i] = np.array([0.0,0.0,0.0])\n",
    "            insula_output_terminal[i] = np.array([0.0,0.0,0.0])\n",
    "        if cortex_neuron[neuron]['soma_position'][0] < IC_a_m:\n",
    "                insula_output[i][0] += j\n",
    "        else:\n",
    "            if cortex_neuron[neuron]['soma_position'][0] < IC_m_p:\n",
    "                insula_output[i][1] += j\n",
    "            else:\n",
    "                insula_output[i][2] += j\n",
    "    if cortex_neuron[neuron]['soma_position'][0] < IC_a_m:\n",
    "        insula_output_count[0] += 1\n",
    "    else:\n",
    "        if cortex_neuron[neuron]['soma_position'][0] < IC_m_p:\n",
    "            insula_output_count[1] += 1\n",
    "        else:\n",
    "            insula_output_count[2] += 1\n",
    "for neuron in cortex_neuron:\n",
    "    for i,j in zip(target_region,cortex_neuron[neuron]['axonal_projection_terminal_combined']):\n",
    "        if cortex_neuron[neuron]['soma_position'][0] < IC_a_m:\n",
    "                insula_output_terminal[i][0] += j\n",
    "        else:\n",
    "            if cortex_neuron[neuron]['soma_position'][0] < IC_m_p:\n",
    "                insula_output_terminal[i][1] += j\n",
    "            else:\n",
    "                insula_output_terminal[i][2] += j\n",
    "for i in insula_output:\n",
    "    insula_output[i] /= insula_output_count\n",
    "    insula_output_terminal[i] /= insula_output_count"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.save(data_path+'insula_input.npy',insula_input)\n",
    "np.save(data_path+'insula_input_terminal.npy',insula_input_terminal)\n",
    "np.save(data_path+'insula_output.npy',insula_output)\n",
    "np.save(data_path+'insula_output_terminal.npy',insula_output_terminal)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#import networkx as nx\n",
    "#from netgraph import Graph\n",
    "node_layout = {\n",
    "    'ORB':np.array([-4.25-1,0.25]),'ILA':np.array([-3.5-1,1.75]),'PL':np.array([-3.5-1,1.25]),'FRP':np.array([-3.75-1,1]),'ACA':np.array([-3.0-1,1.75]),\n",
    "    'MOs':np.array([-2-1,2]),'MOp':np.array([-1.5-1,2]),'SSp':np.array([-1.0-1,2]),'SSs':np.array([-0.5-1,2]),\n",
    "    'TEa':np.array([-0.0-1,2]),'ECT':np.array([0.5-1,2]),'PERI':np.array([1-1,2]),\n",
    "    'ICa':np.array([-0.4,-0.2]),'ICm':np.array([0,0.45]),'ICp':np.array([0.4,-0.2]),\n",
    "    'PIR':np.array([-3,-2]),'TR':np.array([-2.5,-2.0]),'CTXsp':np.array([-2,-2]),'RHP':np.array([-1,-1.5]),\n",
    "    'PAL':np.array([0.25,-2]),'STRd':np.array([1.5,-1.5]),'STRv':np.array([1,-2]),'sAMY':np.array([-0.5,-2]),'COA':np.array([-1,-2]),\n",
    "    'MED':np.array([2,2]),'ILM':np.array([2,1.5]),'VENT':np.array([2,1]),\n",
    "    'LZ':np.array([3,-0.2]),\n",
    "    'PAG':np.array([4,-0.5]),'MRN':np.array([4,-1]),\n",
    "    'PB':np.array([6,0.5]),\n",
    "    'PARN':np.array([6,0]),'NTS':np.array([6.5,0.25]),\n",
    "}\n",
    "threshold = 500"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## IC input"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mode = 0\n",
    "color_mode = 1\n",
    "insula_network = nx.DiGraph()\n",
    "node_name = ['ICa','ICm','ICp']\n",
    "node_labels = {'ICa':'ICa','ICm':'ICm','ICp':'ICp'}\n",
    "node_size = {'ICa':15,'ICm':15,'ICp':15}\n",
    "node_color = {'ICa':'#DB4537','ICm':'#55AF7B','ICp':'#4286F3'}\n",
    "node_color_count = {}\n",
    "edge_color = {}\n",
    "edge_width = {}\n",
    "node_sets = {'input':list(tmp_transform_1.keys())+['ICa','ICm','ICp'],'output':target_region}\n",
    "for i in node_name:\n",
    "    insula_network.add_node(i)\n",
    "for i in insula_input:\n",
    "    if i in ['ICa','ICm','ICp']:\n",
    "        continue\n",
    "    tmp = 0\n",
    "    for j,k in zip(['ICa','ICm','ICp'],insula_input[i]):\n",
    "        if k<threshold:\n",
    "            continue\n",
    "        if i not in node_name:\n",
    "            node_name.append(i)\n",
    "            insula_network.add_node(i)\n",
    "            if color_mode==0:\n",
    "                node_color[i] = '#fac230'\n",
    "            else:\n",
    "                node_color[i] = np.array([0.0,0.0,0.0])\n",
    "        insula_network.add_edge(i,j)\n",
    "        if mode==0:\n",
    "            edge_width[i,j] = k/1000\n",
    "        else:\n",
    "            edge_width[i,j] = np.log(k)/3\n",
    "        tmp += edge_width[i,j]\n",
    "        edge_color[i,j] = node_color[j]\n",
    "        if color_mode==1:\n",
    "            node_color[i][['ICa','ICm','ICp'].index(j)] += edge_width[i,j]\n",
    "    node_size[i] = tmp**0.5*5\n",
    "if color_mode==1:\n",
    "    for i in node_color:\n",
    "        if i in ['ICa','ICm','ICp']:\n",
    "            continue\n",
    "        factors = node_color[i]/np.sum(node_color[i])\n",
    "        node_color[i] = factors[0]*np.array([0.859,0.267,0.216])+factors[1]*np.array([0.059,0.616,0.345])+factors[2]*np.array([0.258,0.52,0.957])\n",
    "        node_color[i] = rgb_to_hex(np.round(node_color[i]*255).astype(int))\n",
    "for i in insula_input:\n",
    "    if i not in ['ICa','ICm','ICp']:\n",
    "        continue\n",
    "    for j,k in zip(['ICa','ICm','ICp'],insula_input[i]):\n",
    "        if k<threshold:\n",
    "            continue\n",
    "        insula_network.add_edge(i,j)\n",
    "        if mode==0:\n",
    "            edge_width[i,j] = k/1000\n",
    "        else:\n",
    "            edge_width[i,j] = np.log(k)/3\n",
    "        edge_color[i,j] = node_color[i]\n",
    "fig = plt.figure(figsize=(12,12))\n",
    "ax = plt.subplot(1,1,1)\n",
    "tmp = Graph(insula_network,\n",
    "      node_size=node_size,\n",
    "      node_color=node_color, \n",
    "      node_edge_width=1,\n",
    "      node_labels=True,\n",
    "      node_label_fontdict=dict(size=7),\n",
    "      node_label_offset=0.0,\n",
    "      edge_color=edge_color,  \n",
    "      edge_width=edge_width,        \n",
    "      edge_alpha=0.75,\n",
    "      node_layout = node_layout,\n",
    "      edge_layout='bundled',\n",
    "      edge_layout_kwargs={'k':250,'total_cycles':8,'total_iterations':50,'compatibility_threshold':0.0025},\n",
    "      with_labels=True,\n",
    "      ax=ax,\n",
    ")\n",
    "ax.set_xlim(-6,7)\n",
    "ax.set_ylim(-3,3)\n",
    "plt.savefig(data_path+'IC_input.svg',format='svg')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## IC output"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mode = 0\n",
    "color_mode = 1\n",
    "insula_network = nx.DiGraph()\n",
    "node_name = ['ICa','ICm','ICp']\n",
    "node_labels = {'ICa':'ICa','ICm':'ICm','ICp':'ICp'}\n",
    "node_size = {'ICa':15,'ICm':15,'ICp':15}\n",
    "node_color = {'ICa':'#DB4537','ICm':'#55AF7B','ICp':'#4286F3'}\n",
    "edge_color = {}\n",
    "edge_width = {}\n",
    "node_sets = {'input':list(tmp_transform_1.keys())+['ICa','ICm','ICp'],'output':target_region}\n",
    "for i in node_name:\n",
    "    insula_network.add_node(i)\n",
    "for i in insula_input:\n",
    "    if i not in ['ICa','ICm','ICp']:\n",
    "        continue\n",
    "    for j,k in zip(['ICa','ICm','ICp'],insula_input[i]):\n",
    "        if k<500:\n",
    "            continue\n",
    "        insula_network.add_edge(i,j)\n",
    "        if mode==0:\n",
    "            edge_width[i,j] = k/1000\n",
    "        else:\n",
    "            edge_width[i,j] = np.log(k)/3\n",
    "        edge_color[i,j] = node_color[i]\n",
    "for i in insula_output:\n",
    "    if i in ['AI','GU','VISC']:\n",
    "        continue\n",
    "    tmp = 0\n",
    "    for j,k in zip(['ICa','ICm','ICp'],insula_output[i]):\n",
    "        if k<500:\n",
    "            continue\n",
    "        if i not in node_name:\n",
    "            node_name.append(i)\n",
    "            if color_mode==0:\n",
    "                node_color[i] = '#f032e6'\n",
    "            else:\n",
    "                node_color[i] = np.array([0.0,0.0,0.0])\n",
    "            insula_network.add_node(i)\n",
    "        insula_network.add_edge(j,i)\n",
    "        if mode==0:\n",
    "            edge_width[j,i] = k/1000\n",
    "            if i=='STRd':\n",
    "                edge_width[j,i] /= 1\n",
    "        else:\n",
    "            edge_width[j,i] = np.log(k)/3\n",
    "        tmp += edge_width[j,i]\n",
    "        edge_color[j,i] = node_color[j]\n",
    "        if color_mode==1:\n",
    "            node_color[i][['ICa','ICm','ICp'].index(j)] += edge_width[j,i]\n",
    "    node_size[i] = tmp**0.5*5\n",
    "if color_mode==1:\n",
    "    for i in node_color:\n",
    "        if i in ['ICa','ICm','ICp']:\n",
    "            continue\n",
    "        factors = node_color[i]/np.sum(node_color[i])\n",
    "        node_color[i] = factors[0]*np.array([0.859,0.267,0.216])+factors[1]*np.array([0.059,0.616,0.345])+factors[2]*np.array([0.258,0.52,0.957])\n",
    "        node_color[i] = rgb_to_hex(np.round(node_color[i]*255).astype(int))\n",
    "fig = plt.figure(figsize=(12,12))\n",
    "ax = plt.subplot(1,1,1)\n",
    "tmp = Graph(insula_network,\n",
    "      node_size=node_size,\n",
    "      node_color=node_color, \n",
    "      node_edge_width=1,\n",
    "      node_labels=True,\n",
    "      node_label_fontdict=dict(size=7),\n",
    "      node_label_offset=0.0,\n",
    "      edge_color=edge_color,  \n",
    "      edge_width=edge_width,        \n",
    "      edge_alpha=0.75,        \n",
    "      node_layout = node_layout,\n",
    "      edge_layout='bundled',\n",
    "      edge_layout_kwargs={'k':250,'total_cycles':8,'total_iterations':50,'compatibility_threshold':0.0025},\n",
    "      with_labels=True,\n",
    "      ax=ax,\n",
    ")\n",
    "ax.set_xlim(-6,7)\n",
    "ax.set_ylim(-3,3)\n",
    "plt.savefig(data_path+'IC_ouput.svg',format='svg')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## IC in and output"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mode = 0\n",
    "color_mode = 1\n",
    "insula_network = nx.DiGraph()\n",
    "node_name = ['ICa','ICm','ICp']\n",
    "node_labels = {'ICa':'ICa','ICm':'ICm','ICp':'ICp'}\n",
    "node_size = {'ICa':15,'ICm':15,'ICp':15}\n",
    "node_color = {'ICa':'#DB4537','ICm':'#55AF7B','ICp':'#4286F3'}\n",
    "edge_color = {}\n",
    "edge_width = {}\n",
    "node_sets = {'input':list(tmp_transform_1.keys())+['ICa','ICm','ICp'],'output':target_region}\n",
    "for i in node_name:\n",
    "    insula_network.add_node(i)\n",
    "input_node_tmp = []\n",
    "for i in insula_input:\n",
    "    if i in ['ICa','ICm','ICp']:\n",
    "        continue\n",
    "    tmp = 0\n",
    "    mark = 0\n",
    "    for j,k in zip(['ICa','ICm','ICp'],insula_input[i]):\n",
    "        if k<500:\n",
    "            continue\n",
    "        if i not in node_name:\n",
    "            node_name.append(i)\n",
    "            insula_network.add_node(i)\n",
    "            if color_mode==0:\n",
    "                node_color[i] = '#fac230'\n",
    "            else:\n",
    "                node_color[i] = np.array([0.0,0.0,0.0])\n",
    "        insula_network.add_edge(i,j)\n",
    "        if mode==0:\n",
    "            edge_width[i,j] = k/1000\n",
    "        else:\n",
    "            edge_width[i,j] = np.log(k)/3\n",
    "        tmp += edge_width[i,j]\n",
    "        edge_color[i,j] = '#fac230'\n",
    "        if color_mode==1:\n",
    "            node_color[i][['ICa','ICm','ICp'].index(j)] += edge_width[i,j]\n",
    "        mark = 1\n",
    "    if mark:\n",
    "        input_node_tmp.append(i)\n",
    "        node_size[i] = tmp\n",
    "for i in insula_input:\n",
    "    if i not in ['ICa','ICm','ICp']:\n",
    "        continue\n",
    "    for j,k in zip(['ICa','ICm','ICp'],insula_input[i]):\n",
    "        if k<500:\n",
    "            continue\n",
    "        insula_network.add_edge(i,j)\n",
    "        if mode==0:\n",
    "            edge_width[i,j] = k/1000\n",
    "        else:\n",
    "            edge_width[i,j] = np.log(k)/3\n",
    "        edge_color[i,j] = node_color[i]\n",
    "for i in insula_output:\n",
    "    if i in ['AI','GU','VISC']:\n",
    "        continue\n",
    "    tmp = 0\n",
    "    for j,k in zip(['ICa','ICm','ICp'],insula_output[i]):\n",
    "        if k<500:\n",
    "            continue\n",
    "        if i not in node_name:\n",
    "            node_name.append(i)\n",
    "            insula_network.add_node(i)\n",
    "            if color_mode==0:\n",
    "                node_color[i] = '#f032e6'\n",
    "            else:\n",
    "                if i in input_node_tmp:\n",
    "                    node_color[i] += insula_output[i]\n",
    "                else:\n",
    "                    node_color[i] = np.array([0.0,0.0,0.0])\n",
    "        insula_network.add_edge(j,i)\n",
    "        if mode==0:\n",
    "            edge_width[j,i] = k/1000\n",
    "            if i=='STRd':\n",
    "                edge_width[j,i] /= 1\n",
    "        else:\n",
    "            edge_width[j,i] = np.log(k)/3\n",
    "        tmp += edge_width[j,i]\n",
    "        edge_color[j,i] = node_color[j]\n",
    "        if color_mode==1:\n",
    "            node_color[i][['ICa','ICm','ICp'].index(j)] += edge_width[j,i]\n",
    "    if i in input_node_tmp:\n",
    "        node_size[i] += tmp\n",
    "    else:\n",
    "        node_size[i] = tmp**0.5*5\n",
    "for i in input_node_tmp:\n",
    "    node_size[i] = node_size[i]**0.5*5\n",
    "    factors = node_color[i]/np.sum(node_color[i])\n",
    "if color_mode==1:\n",
    "    for i in node_color:\n",
    "        if i in ['ICa','ICm','ICp']:\n",
    "            continue\n",
    "        factors = node_color[i]/np.sum(node_color[i])\n",
    "        node_color[i] = factors[0]*np.array([0.859,0.267,0.216])+factors[1]*np.array([0.059,0.616,0.345])+factors[2]*np.array([0.258,0.52,0.957])\n",
    "        node_color[i] = rgb_to_hex(np.round(node_color[i]*255).astype(int))\n",
    "        for j in edge_color:\n",
    "            if j[0]==i:\n",
    "                edge_color[j] = node_color[i]\n",
    "fig = plt.figure(figsize=(12,12))\n",
    "ax = plt.subplot(1,1,1)\n",
    "tmp = Graph(insula_network,\n",
    "      node_size=node_size,\n",
    "      node_color=node_color, \n",
    "      node_edge_width=1,\n",
    "      node_labels=True,\n",
    "      node_label_fontdict=dict(size=7),\n",
    "      node_label_offset=0.0,\n",
    "      edge_color=edge_color,  \n",
    "      edge_width=edge_width,        \n",
    "      edge_alpha=0.75,        \n",
    "      #node_layout='spring',\n",
    "      #node_layout='circular',\n",
    "      #node_layout_kwargs=dict(subsets=node_sets),\n",
    "      node_layout = node_layout,\n",
    "      edge_layout='bundled',\n",
    "      edge_layout_kwargs={'k':250,'total_cycles':8,'total_iterations':50,'compatibility_threshold':0.0025},\n",
    "      with_labels=True,\n",
    "      ax=ax,\n",
    ")\n",
    "ax.set_xlim(-6,7)\n",
    "ax.set_ylim(-3,3)\n",
    "plt.savefig(data_path+'IC_inouput.svg',format='svg')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## iso ic networks qutaty"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ISO_IC_projection_distribution = {}\n",
    "for i in target_region_hierarchy['Cortex']+['ICa','ICm','ICp']:\n",
    "    if i in ['AI','GU','VISC']:\n",
    "        continue\n",
    "    if i in ['ICa','ICm','ICp']:\n",
    "        index = ['ICa','ICm','ICp'].index(i)\n",
    "        ISO_IC_projection_distribution[i] = np.array([insula_input[i],insula_input[i]])\n",
    "        ISO_IC_projection_distribution[i][1] = np.array([insula_input[j][index] for j in ['ICa','ICm','ICp']])\n",
    "    else:\n",
    "        ISO_IC_projection_distribution[i] = np.array([insula_input[i],insula_output[i]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "name_index = list(ISO_IC_projection_distribution.keys())\n",
    "_iso_to_ic = np.array([[ISO_IC_projection_distribution[i][0][0] for i in name_index],[ISO_IC_projection_distribution[i][0][1] for i in name_index],[ISO_IC_projection_distribution[i][0][2] for i in name_index]])\n",
    "_ic_to_iso = np.array([[ISO_IC_projection_distribution[i][1][0] for i in name_index],[ISO_IC_projection_distribution[i][1][1] for i in name_index],[ISO_IC_projection_distribution[i][1][2] for i in name_index]])\n",
    "plt.ylim(-40000,40000)\n",
    "if 0:\n",
    "    _iso_to_ic = np.log(_iso_to_ic+1)\n",
    "    _ic_to_iso = np.log(_ic_to_iso+1)\n",
    "    plt.ylim(-30,30)\n",
    "# iso_to_ic\n",
    "tmp = 0\n",
    "for i in range(2,-1,-1):\n",
    "    plt.bar(range(len(name_index)),_iso_to_ic[i]+tmp,color=region_transform_color[['AIa','AIm','AIp'][i]],zorder=i)\n",
    "    tmp += _iso_to_ic[i]\n",
    "# ic_to_iso\n",
    "tmp = 0\n",
    "for i in range(2,-1,-1):\n",
    "    plt.bar(range(len(name_index)),-(_ic_to_iso[i]+tmp),color=region_transform_color[['AIa','AIm','AIp'][i]],zorder=i)\n",
    "    tmp += _ic_to_iso[i]\n",
    "_ = plt.xticks(range(len(name_index)),name_index,rotation=-90)\n",
    "plt.axhline(0,color='k')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tmp = np.sum(_iso_to_ic,axis=0)/(np.sum(_iso_to_ic,axis=0)+np.sum(_ic_to_iso,axis=0))\n",
    "plt.bar(range(len(name_index)),tmp,width=0.5,color='#fa3c32')\n",
    "_ = plt.xticks(range(len(name_index)),name_index,rotation=-90)\n",
    "for i in range(len(name_index)):\n",
    "    plt.text(i,tmp[i],str(np.round(tmp[i],1)), ha='center', va='bottom')\n",
    "plt.ylim(0,1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## NTS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "index = target_region.index('NTS')\n",
    "flatmap = cv.imread('../resource/flatmapedge2.tif',1)\n",
    "fig = plt.figure(figsize=(5,5))\n",
    "ax = fig.add_subplot(1,1,1)\n",
    "for i in insular_neuron_name:\n",
    "    if cortex_neuron[i]['axonal_projection_length_combined'][index]<0:\n",
    "        continue\n",
    "    position = cortex_neuron[i]['soma_flatten_position']\n",
    "    if np.sum(abs(np.array(position)))!=0:\n",
    "        cv.circle(flatmap,[int(position[0])*2, int(position[1])*2],int(np.sqrt(cortex_neuron[i]['axonal_projection_length_combined'][index]))//5,[0,255,0],-1)\n",
    "ax.matshow(flatmap)\n",
    "ax.set_xticks([],[])\n",
    "ax.set_yticks([],[])\n",
    "#ax.set_axis_off()\n",
    "ax.set_xlabel('Lat. to Med. axis',loc='center',fontsize=fontsize_label*0.6)\n",
    "ax.set_ylabel('Ant. to Post. axis',rotation=-90,loc='center',fontsize=fontsize_label*0.6).set_verticalalignment('top')\n",
    "ax.set_title('Soma distribution of the NTS-projection neurons',fontsize=fontsize_title*0.5)\n",
    "ax.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = []\n",
    "for i in insular_neuron_name:\n",
    "    if cortex_neuron[i]['cortex_subregion']=='d':\n",
    "        x.append((5472-cortex_neuron[i]['soma_position'][0])/1000)\n",
    "np.min(x),np.max(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(1.54-1.21)/(3.31-1.21)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = []\n",
    "for i in insular_neuron_name:\n",
    "    if cortex_neuron[i]['cortex_subregion']=='p':\n",
    "        x.append((5472-cortex_neuron[i]['soma_position'][0])/1000)\n",
    "np.min(x),np.max(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "index = target_region.index('NTS')\n",
    "info = [[],[]]\n",
    "fig = plt.figure(figsize=(5,5))\n",
    "ax = fig.add_subplot(1,1,1)\n",
    "for i in insular_neuron_name:\n",
    "    if cortex_neuron[i]['axonal_projection_length_combined'][index]<0:\n",
    "        continue\n",
    "    position = (5472-cortex_neuron[i]['soma_position'][0])/1000\n",
    "    strength = np.sqrt(cortex_neuron[i]['axonal_projection_length_combined'][index])**2\n",
    "    info[0].append(position)\n",
    "    info[1].append(strength)\n",
    "for i,j in zip(info[0],info[1]):\n",
    "    ax.plot([i,i],[0,j],color='#55af7b',linewidth=0.3)\n",
    "#ax.set_xticks([],[])\n",
    "ax.set_yticks([],[])\n",
    "#ax.set_axis_off()\n",
    "ax.set_xlabel('Bregma',loc='center',fontsize=fontsize_label*0.6)\n",
    "ax.set_ylabel('Strength',rotation=-90,loc='center',fontsize=fontsize_label*0.6).set_verticalalignment('top')\n",
    "ax.set_title('Soma distribution of the NTS-projection neurons',fontsize=fontsize_title*0.5)\n",
    "ax.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "selected_region = ['CEA','BLA']\n",
    "insula_extra_projection_info = {}\n",
    "for i in tool_bar(insular_neuron_name):\n",
    "    insula_extra_projection_info[i] = {}\n",
    "    control_console_output(0)\n",
    "    for j,k in zip(selected_region,get_axonal_statics_by_regions(i,selected_region,'projectregion')):\n",
    "        insula_extra_projection_info[i][j] = k\n",
    "    control_console_output(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for r in selected_region:\n",
    "    flatmap = cv.imread('../resource/flatmapedge2.tif',1)\n",
    "    fig = plt.figure(figsize=(5,5))\n",
    "    ax = fig.add_subplot(1,1,1)\n",
    "    for i in insular_neuron_name:\n",
    "        if insula_extra_projection_info[i][r]<0:\n",
    "            continue\n",
    "        position = cortex_neuron[i]['soma_flatten_position']\n",
    "        if np.sum(abs(np.array(position)))!=0:\n",
    "            cv.circle(flatmap,[int(position[0])*2, int(position[1])*2],int(np.sqrt(insula_extra_projection_info[i][r]))//10,[15,157,88],-1)\n",
    "    ax.matshow(flatmap)\n",
    "    ax.set_xticks([],[])\n",
    "    ax.set_yticks([],[])\n",
    "    #ax.set_axis_off()\n",
    "    ax.set_xlabel('Lat. to Med. axis',loc='center',fontsize=fontsize_label*0.6)\n",
    "    ax.set_ylabel('Ant. to Post. axis',rotation=-90,loc='center',fontsize=fontsize_label*0.6).set_verticalalignment('top')\n",
    "    ax.set_title('Soma distribution of the'+r+'-projection neurons',fontsize=fontsize_title*0.5)\n",
    "    ax.legend()\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Dendrite Analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_neuron_tree(neuron):\n",
    "    return iondata.getNeuronTreeByID(neuron[:-7],neuron[-7:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if not os.path.exists(data_path+'insular_dendrite_info.npy'):\n",
    "    insular_dendrite_info = {}\n",
    "    for neuron in tool_bar(insular_neuron_name):\n",
    "        neuron_dendrites = get_neuron_tree(neuron).getDendrite()\n",
    "        if len(neuron_dendrites)==0:\n",
    "            continue\n",
    "        insular_dendrite_info[neuron] = []\n",
    "        for edge in neuron_dendrites:\n",
    "            insular_dendrite_info[neuron].append([get_mirrored_point(point.xyz) for point in edge.data])\n",
    "    np.save(data_path+'insular_dendrite_info.npy',insular_dendrite_info)\n",
    "else:\n",
    "    insular_dendrite_info = np.load(data_path+'insular_dendrite_info.npy',allow_pickle=True).item()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## dendrites info -> features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if not os.path.exists(data_path+'insular_dendrite_feature.npy'):\n",
    "    insular_dendrite_feature = {}\n",
    "    for neuron in tool_bar(insular_dendrite_info):\n",
    "        insular_dendrite_feature[neuron] = {}\n",
    "        insular_dendrite_feature[neuron]['soma_position'] = cortex_neuron[neuron]['soma_flatten_position']\n",
    "        insular_dendrite_feature[neuron]['soma_depth'] = cortex_neuron[neuron]['soma_depth']\n",
    "        insular_dendrite_feature[neuron]['soma_outter_position'] = cortex_neuron[neuron]['soma_outter_position']\n",
    "        insular_dendrite_feature[neuron]['soma_inner_position'] = cortex_neuron[neuron]['soma_inner_position']\n",
    "        insular_dendrite_feature[neuron]['dendrite_position'] = []\n",
    "        insular_dendrite_feature[neuron]['dendrite_depth'] = []\n",
    "        insular_dendrite_feature[neuron]['dendrite_outter_position'] = []\n",
    "        insular_dendrite_feature[neuron]['dendrite_inner_position'] = []\n",
    "        for edge in insular_dendrite_info[neuron]:\n",
    "            edge_position = []\n",
    "            edge_depth = []\n",
    "            edge_outter_position = []\n",
    "            edge_inner_position = []\n",
    "            for i in compute_flatten_position_and_depth(edge[::20]+[edge[-1]]):\n",
    "                edge_depth.append(i[0])\n",
    "                edge_outter_position.append(i[1])\n",
    "                edge_inner_position.append(i[2])\n",
    "                edge_position.append(i[3])\n",
    "            insular_dendrite_feature[neuron]['dendrite_position'].append(edge_position)\n",
    "            insular_dendrite_feature[neuron]['dendrite_depth'].append(edge_depth)\n",
    "            insular_dendrite_feature[neuron]['dendrite_outter_position'].append(edge_outter_position)\n",
    "            insular_dendrite_feature[neuron]['dendrite_inner_position'].append(edge_inner_position)\n",
    "    np.save(data_path+'insular_dendrite_feature.npy',insular_dendrite_feature)\n",
    "else:\n",
    "    insular_dendrite_feature = np.load(data_path+'insular_dendrite_feature.npy',allow_pickle=True).item()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "kernel_sigma = 3\n",
    "_min,_max = 0,1200\n",
    "bins = 100\n",
    "step = (_max-_min)/(bins+1)\n",
    "distribution_index = np.array([i for i in range(bins)])\n",
    "neuron_name_dendrite = []\n",
    "feature_matrix = []\n",
    "for i in insular_dendrite_feature:\n",
    "    if insular_dendrite_feature[i]['soma_depth']==0:\n",
    "        continue\n",
    "    neuron_name_dendrite.append(i)\n",
    "    feature_matrix.append(np.zeros(bins,))\n",
    "    for j in insular_dendrite_feature[i]['dendrite_depth']:\n",
    "        for k in j:\n",
    "            if k==0:\n",
    "                continue\n",
    "            center = k*20/step\n",
    "            feature_matrix[-1] += func_gaussian(distribution_index,center=center,sigma=kernel_sigma)\n",
    "    if 0:\n",
    "        feature_matrix[-1] /= np.sum(feature_matrix[-1])\n",
    "feature_matrix = np.array(feature_matrix)\n",
    "feature_matrix = np.log(feature_matrix+1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "kernel_sigma = 5\n",
    "_min,_max = 0,1200\n",
    "bins = 100\n",
    "r = 10\n",
    "r2 = 5\n",
    "step = (_max-_min)/(bins+1)\n",
    "distribution_index = np.array([i for i in range(bins)])\n",
    "neuron_name_dendrite = []\n",
    "feature_matrix = []\n",
    "processed_feature_matrix = []\n",
    "for i in insular_dendrite_feature:\n",
    "    if insular_dendrite_feature[i]['soma_depth']==0:\n",
    "        continue\n",
    "    neuron_name_dendrite.append(i)\n",
    "    feature_matrix.append(np.zeros(bins,))\n",
    "    processed_feature_matrix.append(np.zeros(5,))\n",
    "    for j in insular_dendrite_feature[i]['dendrite_depth']:\n",
    "        for k in j:\n",
    "            if k==0:\n",
    "                continue\n",
    "            center = k*20/step\n",
    "            feature_matrix[-1] += func_gaussian(distribution_index,center=center,sigma=kernel_sigma)\n",
    "    processed_feature_matrix[-1][0] = np.sum(feature_matrix[-1][:r]) #L1\n",
    "    soma_l = int(np.round(insular_dendrite_feature[i]['soma_depth']*20/step))\n",
    "    processed_feature_matrix[-1][1] = np.sum(feature_matrix[-1][np.max([soma_l-2*r2,r]):np.max([soma_l-r2,r])]) #soma up2\n",
    "    processed_feature_matrix[-1][2] = np.sum(feature_matrix[-1][np.max([soma_l-r2,r]):np.max([soma_l,r])]) #soma up\n",
    "    processed_feature_matrix[-1][3] = np.sum(feature_matrix[-1][np.max([soma_l,r]):np.max([soma_l+r2,r])]) #soma down\n",
    "    processed_feature_matrix[-1][4] = np.sum(feature_matrix[-1][np.max([soma_l+r2,r]):np.max([soma_l+2*r2,r])]) #soma down2\n",
    "feature_matrix = np.array(processed_feature_matrix)\n",
    "feature_matrix = np.log(feature_matrix+1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#overall information of the clustering\n",
    "if 0:\n",
    "    tsne = TSNE(n_components=2, perplexity=35, random_state=0,n_iter=10000)\n",
    "    X_tsne = tsne.fit_transform(feature_matrix)\n",
    "hierarchy.set_link_color_palette(cluster_color_hex)\n",
    "color_threhold = 30\n",
    "fig = plt.figure(figsize=(10,5))\n",
    "if not os.path.exists(cluster_path+'dendrite_cluster_object2.npy'):\n",
    "    cluster_object = hierarchy.dendrogram(hierarchy.linkage(feature_matrix,'ward'),color_threshold=color_threhold,distance_sort=True)\n",
    "    np.save(cluster_path+'cluster_object.npy',cluster_object)\n",
    "else:\n",
    "    cluster_object = np.load(cluster_path+'dendrite_cluster_object.npy',allow_pickle=True).item()\n",
    "plt.axhline(color_threhold, linestyle='dashed', color='#AAAAAA')\n",
    "plt.title(\"Hierachy of the clustering\",fontsize=20)\n",
    "plt.locator_params(axis='y', nbins=3)\n",
    "plt.xticks([],[])\n",
    "plt.xlabel('neurons = '+str(len(neuron_name_dendrite)),fontsize=15)\n",
    "plt.ylabel('Cluster distance',fontsize=15)\n",
    "#plt.savefig(cluster_path+'cluster_hierachy.svg',format='svg')\n",
    "plt.show()\n",
    "cluster_object_parse = [[],[]]\n",
    "for i,j in zip(cluster_object['leaves'],cluster_object['leaves_color_list']):\n",
    "    neuron = neuron_name_dendrite[i]\n",
    "    projection_matrix = feature_matrix[i]\n",
    "    cluster_object_parse[0].append(cluster_color_hex.index(j))\n",
    "    cluster_object_parse[1].append(projection_matrix)\n",
    "cluster_object_parse[1] = np.array(cluster_object_parse[1])\n",
    "cluster_object_parse_visualization = np.zeros(shape=(len(neuron_name_dendrite),bins,3))\n",
    "cluster_object_parse_visualization[:,0,:] = np.array([cluster_color[i] for i in cluster_object_parse[0]])\n",
    "tmp = np.array([cluster_object_parse[1],cluster_object_parse[1]*0,cluster_object_parse[1]*0]).transpose((1,2,0))\n",
    "tmp /= np.max(tmp)\n",
    "cluster_object_parse_visualization = tmp*1\n",
    "fig = plt.figure(figsize=(4,6))\n",
    "ax = fig.add_subplot()\n",
    "ax.matshow(cluster_object_parse_visualization.transpose((1,0,2)),vmax=0.05)\n",
    "_ = ax.set_xticks([],[])\n",
    "ax.set_xlabel('neurons = '+str(len(neuron_name_dendrite)),fontsize=fontsize_label)\n",
    "_ = ax.set_yticks([],[])\n",
    "#plt.savefig(cluster_path+'cluster_projection_matrix.svg',format='svg')\n",
    "plt.show()\n",
    "tmp = {}\n",
    "for i,j in enumerate(neuron_name_dendrite):\n",
    "    _cluster = cluster_color_hex.index(cluster_object['leaves_color_list'][cluster_object['leaves'].index(i)])\n",
    "    if _cluster not in tmp:\n",
    "        tmp[_cluster] = []\n",
    "    tmp[_cluster].append([j,cortex_neuron[j]['soma_flatten_position'],X_tsne[i]])\n",
    "_flatmap = cv.imread('../resource/flatmapedge2.tif',1)\n",
    "fig = plt.figure(figsize=(10,10))\n",
    "flatmap = _flatmap*1\n",
    "for i in range(len(list(tmp.keys()))):\n",
    "    tmp2 = []\n",
    "    tmp3 = []\n",
    "    for j in tmp[i]:\n",
    "        if np.sum(abs(np.array(j[1])))!=0:\n",
    "            cv.circle(flatmap,[int(j[1][0])*2, int(j[1][1])*2],7,cluster_color[i]*255,-1)\n",
    "        tmp2.append(j[2])\n",
    "        tmp3.append(j[0])\n",
    "    tmp2 = np.array(tmp2)\n",
    "    plt.scatter(tmp2[:,0], tmp2[:,1],c=cluster_color_hex[i],s=30,marker='o',label=str(i+1))\n",
    "    save_scene(tmp3,cluster_color[i],cluster_path+'dendrite_cluster_'+str(i+1))\n",
    "plt.legend()\n",
    "plt.title(\"t-SNE Visualization of the clusters\",fontsize=fontsize_title)\n",
    "plt.axis('off')\n",
    "#plt.savefig(cluster_path+'cluster_t-SNE.svg',format='svg')\n",
    "plt.show()\n",
    "fig = plt.figure(figsize=(10,10))\n",
    "plt.imshow(flatmap)\n",
    "plt.xlim(300,2700)\n",
    "plt.ylim(2700,300)\n",
    "plt.title(\"The soma distribution of the clusters\",fontsize=fontsize_title)\n",
    "plt.axis('off')\n",
    "#plt.savefig(cluster_path+'cluster_distribution.svg',format='svg')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "random_sample = 3\n",
    "sample_x_step = 300\n",
    "for i in range(5):\n",
    "    fig = plt.figure(figsize=(10,4))\n",
    "    bias = 0\n",
    "    q = []\n",
    "    for j in np.random.choice(range(len(tmp[i])),random_sample):\n",
    "        neuron = tmp[i][j][0]\n",
    "        print(neuron)\n",
    "        q.append(feature_matrix[neuron_name_dendrite.index(neuron)])\n",
    "        for k,l in zip(insular_dendrite_feature[neuron]['dendrite_outter_position'],insular_dendrite_feature[neuron]['dendrite_depth']):\n",
    "            x = np.array(k)[:,0]-insular_dendrite_feature[neuron]['soma_outter_position'][0]\n",
    "            y = np.array(l)\n",
    "            mask = (y!=0).astype(bool)\n",
    "            x = x[mask]*10+bias\n",
    "            y = y[mask]*20\n",
    "            plt.plot(x,y,linewidth=0.3,color=cluster_color_hex[i])\n",
    "        bias += sample_x_step\n",
    "    for j in [135.96993306, 457.4311794 , 840.16536782]:\n",
    "        plt.axhline(j,linestyle='dashed',color='#AAAAAA',linewidth=0.5)\n",
    "    #plt.xlim(300,2700)\n",
    "    plt.ylim(1200,0)\n",
    "    plt.axis('off')\n",
    "    plt.title(\"The dendrite of the cluster \"+str(i+1),fontsize=fontsize_title)\n",
    "    #plt.savefig(cluster_path+'cluster_distribution.svg',format='svg')\n",
    "    plt.show()\n",
    "    plt.matshow(np.array(q))\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# neuronnetwork method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = []\n",
    "Y = []\n",
    "for i in neuron_name_dendrite:\n",
    "    Y.append(i)\n",
    "    X.append(np.zeros(shape=(128,208*3,1)))\n",
    "    for j,k in zip(insular_dendrite_feature[i]['dendrite_outter_position'],insular_dendrite_feature[i]['dendrite_depth']):\n",
    "        x = np.array(j)[:,0]-insular_dendrite_feature[i]['soma_outter_position'][0]\n",
    "        y = np.array(k)\n",
    "        mask = (y!=0).astype(bool)\n",
    "        x = (x[mask]*10/2+64).astype(int)\n",
    "        y = (y[mask]*20/2).astype(int)\n",
    "        for m in range(len(x)-1):\n",
    "            cv2.line(X[-1],[y[m],x[m]],[y[m+1],x[m+1]],1,1)\n",
    "    X[-1] = np.where(X[-1]>0,1,0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(X[-3])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(X[-3])\n",
    "plt.show()\n",
    "plt.imshow(autoencoder.predict(x_train[-3].reshape(1,128,624,1))[0])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(X[-2])\n",
    "plt.show()\n",
    "plt.imshow(autoencoder.predict(x_train[-2].reshape(1,128,624,1))[0])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import tensorflow as tf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "input_img = tf.keras.Input(shape=(128, 208*3, 1))\n",
    "\n",
    "x = tf.keras.layers.Conv2D(16, (3, 3), activation='relu', padding='same')(input_img)\n",
    "x = tf.keras.layers.MaxPooling2D((2, 2), padding='same')(x) #64,104*3\n",
    "x = tf.keras.layers.Conv2D(8, (3, 3), activation='relu', padding='same')(x)\n",
    "x = tf.keras.layers.MaxPooling2D((2, 2), padding='same')(x) #32,52*3\n",
    "x = tf.keras.layers.Conv2D(8, (3, 3), activation='relu', padding='same')(x)\n",
    "x = tf.keras.layers.MaxPooling2D((2, 2), padding='same')(x) #16,26*3\n",
    "x = tf.keras.layers.Conv2D(8, (3, 3), activation='relu', padding='same')(x)\n",
    "encoded = tf.keras.layers.MaxPooling2D((2, 2), padding='same')(x) #8,13*3\n",
    "\n",
    "# at this point the representation is (4, 4, 8) i.e. 128-dimensional\n",
    "\n",
    "x = tf.keras.layers.Conv2D(8, (3, 3), activation='relu', padding='same')(encoded)\n",
    "x = tf.keras.layers.UpSampling2D((2, 2))(x)\n",
    "x = tf.keras.layers.Conv2D(8, (3, 3), activation='relu', padding='same')(x)\n",
    "x = tf.keras.layers.UpSampling2D((2, 2))(x)\n",
    "x = tf.keras.layers.Conv2D(8, (3, 3), activation='relu', padding='same')(x)\n",
    "x = tf.keras.layers.UpSampling2D((2, 2))(x)\n",
    "x = tf.keras.layers.Conv2D(16, (3, 3), activation='relu', padding='same')(x)\n",
    "x = tf.keras.layers.UpSampling2D((2, 2))(x)\n",
    "decoded = tf.keras.layers.Conv2D(1, (3, 3), activation='sigmoid', padding='same')(x)\n",
    "\n",
    "autoencoder = tf.keras.Model(input_img, decoded)\n",
    "autoencoder.compile(optimizer='adam', loss='binary_crossentropy')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_train = np.array(X).astype(float)\n",
    "autoencoder.fit(x_train,x_train,epochs=50,batch_size=32,shuffle=True)\n",
    "x = x_train*1\n",
    "for i in range(9):\n",
    "    x = autoencoder.layers[i](x)\n",
    "feature_matrix = x.numpy().reshape(982,-1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#overall information of the clustering\n",
    "if 0:\n",
    "    tsne = TSNE(n_components=2, perplexity=35, random_state=0,n_iter=10000)\n",
    "    X_tsne = tsne.fit_transform(feature_matrix)\n",
    "hierarchy.set_link_color_palette(cluster_color_hex)\n",
    "color_threhold = 200\n",
    "fig = plt.figure(figsize=(10,5))\n",
    "if not os.path.exists(cluster_path+'dendrite_cluster_object2.npy'):\n",
    "    cluster_object = hierarchy.dendrogram(hierarchy.linkage(X_tsne,'ward'),color_threshold=color_threhold,distance_sort=True)\n",
    "    np.save(cluster_path+'cluster_object.npy',cluster_object)\n",
    "else:\n",
    "    cluster_object = np.load(cluster_path+'dendrite_cluster_object.npy',allow_pickle=True).item()\n",
    "plt.axhline(color_threhold, linestyle='dashed', color='#AAAAAA')\n",
    "plt.title(\"Hierachy of the clustering\",fontsize=20)\n",
    "plt.locator_params(axis='y', nbins=3)\n",
    "plt.xticks([],[])\n",
    "plt.xlabel('neurons = '+str(len(neuron_name_dendrite)),fontsize=15)\n",
    "plt.ylabel('Cluster distance',fontsize=15)\n",
    "#plt.savefig(cluster_path+'cluster_hierachy.svg',format='svg')\n",
    "plt.show()\n",
    "cluster_object_parse = [[],[]]\n",
    "for i,j in zip(cluster_object['leaves'],cluster_object['leaves_color_list']):\n",
    "    neuron = neuron_name_dendrite[i]\n",
    "    projection_matrix = feature_matrix[i]\n",
    "    cluster_object_parse[0].append(cluster_color_hex.index(j))\n",
    "    cluster_object_parse[1].append(projection_matrix)\n",
    "cluster_object_parse[1] = np.array(cluster_object_parse[1])\n",
    "cluster_object_parse_visualization = np.zeros(shape=(len(neuron_name_dendrite),bins,3))\n",
    "cluster_object_parse_visualization[:,0,:] = np.array([cluster_color[i] for i in cluster_object_parse[0]])\n",
    "tmp = np.array([cluster_object_parse[1],cluster_object_parse[1]*0,cluster_object_parse[1]*0]).transpose((1,2,0))\n",
    "tmp /= np.max(tmp)\n",
    "cluster_object_parse_visualization = tmp*1\n",
    "fig = plt.figure(figsize=(4,6))\n",
    "ax = fig.add_subplot()\n",
    "ax.matshow(cluster_object_parse_visualization.transpose((1,0,2)),vmax=0.05)\n",
    "_ = ax.set_xticks([],[])\n",
    "ax.set_xlabel('neurons = '+str(len(neuron_name_dendrite)),fontsize=fontsize_label)\n",
    "_ = ax.set_yticks([],[])\n",
    "#plt.savefig(cluster_path+'cluster_projection_matrix.svg',format='svg')\n",
    "plt.show()\n",
    "tmp = {}\n",
    "for i,j in enumerate(neuron_name_dendrite):\n",
    "    _cluster = cluster_color_hex.index(cluster_object['leaves_color_list'][cluster_object['leaves'].index(i)])\n",
    "    if _cluster not in tmp:\n",
    "        tmp[_cluster] = []\n",
    "    tmp[_cluster].append([j,cortex_neuron[j]['soma_flatten_position'],X_tsne[i]])\n",
    "_flatmap = cv.imread('../resource/flatmapedge2.tif',1)\n",
    "fig = plt.figure(figsize=(10,10))\n",
    "flatmap = _flatmap*1\n",
    "for i in range(len(list(tmp.keys()))):\n",
    "    tmp2 = []\n",
    "    tmp3 = []\n",
    "    for j in tmp[i]:\n",
    "        if np.sum(abs(np.array(j[1])))!=0:\n",
    "            cv.circle(flatmap,[int(j[1][0])*2, int(j[1][1])*2],7,cluster_color[i]*255,-1)\n",
    "        tmp2.append(j[2])\n",
    "        tmp3.append(j[0])\n",
    "    tmp2 = np.array(tmp2)\n",
    "    plt.scatter(tmp2[:,0], tmp2[:,1],c=cluster_color_hex[i],s=30,marker='o',label=str(i+1))\n",
    "    save_scene(tmp3,cluster_color[i],cluster_path+'dendrite_cluster_'+str(i+1))\n",
    "plt.legend()\n",
    "plt.title(\"t-SNE Visualization of the clusters\",fontsize=fontsize_title)\n",
    "plt.axis('off')\n",
    "#plt.savefig(cluster_path+'cluster_t-SNE.svg',format='svg')\n",
    "plt.show()\n",
    "fig = plt.figure(figsize=(10,10))\n",
    "plt.imshow(flatmap)\n",
    "plt.xlim(300,2700)\n",
    "plt.ylim(2700,300)\n",
    "plt.title(\"The soma distribution of the clusters\",fontsize=fontsize_title)\n",
    "plt.axis('off')\n",
    "#plt.savefig(cluster_path+'cluster_distribution.svg',format='svg')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "random_sample = 3\n",
    "sample_x_step = 300\n",
    "for i in range(9):\n",
    "    fig = plt.figure(figsize=(10,4))\n",
    "    bias = 0\n",
    "    q = []\n",
    "    for j in np.random.choice(range(len(tmp[i])),random_sample,replace=False):\n",
    "        neuron = tmp[i][j][0]\n",
    "        print(neuron)\n",
    "        q.append(feature_matrix[neuron_name_dendrite.index(neuron)])\n",
    "        for k,l in zip(insular_dendrite_feature[neuron]['dendrite_outter_position'],insular_dendrite_feature[neuron]['dendrite_depth']):\n",
    "            x = np.array(k)[:,0]-insular_dendrite_feature[neuron]['soma_outter_position'][0]\n",
    "            y = np.array(l)\n",
    "            mask = (y!=0).astype(bool)\n",
    "            x = x[mask]*10+bias\n",
    "            y = y[mask]*20\n",
    "            plt.plot(x,y,linewidth=0.3,color=cluster_color_hex[i])\n",
    "        bias += sample_x_step\n",
    "    for j in [135.96993306, 457.4311794 , 840.16536782]:\n",
    "        plt.axhline(j,linestyle='dashed',color='#AAAAAA',linewidth=0.5)\n",
    "    #plt.xlim(300,2700)\n",
    "    plt.ylim(1200,0)\n",
    "    plt.axis('off')\n",
    "    plt.title(\"The dendrite of the cluster \"+str(i+1),fontsize=fontsize_title)\n",
    "    #plt.savefig(cluster_path+'cluster_distribution.svg',format='svg')\n",
    "    plt.show()\n",
    "    plt.matshow(np.array(q))\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.8.13 ('CX')",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.13"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "6684bfb8f990107dc2e41da7bb6829210ca36b987584a9d9d8b153337675e487"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
